Commit 6b880592 authored by mcariou's avatar mcariou
Browse files

update report

parent 624ea0de
......@@ -48,22 +48,31 @@
\section{Introduction}
For detailed summary of the analysis, see README.md file.
For detailed summary of the analysis, see \textit{README.md} file.
~
Briefly, illumina paired-end reads from 2 Ebola strains, Ebola-800 and Ebola-911 were analysed.
The aim of this study is to identify all differences between the sequences of 800 and 911 compared to the reference sequences.
For this, 2 approaches are possible, first, to assemble genomes de novo and the compare genome sequences to reference genome sequence, second, to align reads against the reference genomes and use software which identify variant directly from the read alignement.
After quality check and read cleaning by Trimmomatic, I 1/ assembled both genome de novo using SPADES 2/ mapped read against a reference genome (Bowtie) and called variants using 2 softwares (bcftools and Freebayes) 3/ Reassembled reads using only mapped reads (SPADES).
~
The second assembly using only reads that mapped against the reference genome was of higher quality. I will use this one in this analysis.
~
The aim of this report is to identify and compare the variations found using these different approaches.
\section{Comparison of genome assemblies with reference genome}
Sequences from the reference genome and the 2 SPADES assembly (2nd version, mapped reads) were manually aligned. See report on the GP gene Ebola\_assembly.pdf.
~
\begin{itemize}
\item Load R packages
<<library, results="hide">>=
library(adegenet)
library(pegas)
......@@ -72,45 +81,62 @@ library(mmod)
library(PopGenome)
library(seqinr)
@
Read the alignment.
<<setup, include=FALSE, cache=FALSE, tidy=TRUE>>=
options(tidy=TRUE, width=70)
@
\item Read the alignment.
<<>>=
workdir<-"~/Documents/CIRI_BIBS_projects/2020_09_Reynard/seq_manual_aln/"
file<-paste0(workdir, "Ebola_aln_assemblyorigin3seq.fst")
fasta<-read.fasta(file)
@
\item Pairwise comparison: Which are the variable positions
<<>>=
workdir<-"~/Documents/CIRI_BIBS_projects/2020_09_Reynard/seq_manual_aln/Ebola_aln_assemblyorigin3seq.fst"
#fasta<-read.dna(workdir, format = "fasta", as.matrix = TRUE)
fasta<-read.fasta(workdir)
# variable positions between ref and 800
pos1<-which(as.vector(fasta$AF086833.2)!=as.vector(fasta$`Ebola-800_assemblypostmap_NODE_1_length_18872_cov_119.776_ID_37`))
#a few base missing at the start and at the end, they are not significant cut between pos 38 and 18910
pos1<-which(as.vector(fasta$AF086833.2)!=
as.vector(fasta$`Ebola-800_assemblypostmap_NODE_1_length_18872_cov_119.776_ID_37`))
# a few base missing at the start and at the end,
# they are not significant cut between pos 38 and 18910
pos1<-pos1[pos1>38 & pos1<18910]
# variable positions between ref and 911
pos2<-which(as.vector(fasta$AF086833.2)!=as.vector(fasta$`C_Ebola-911_assemblypostmap_NODE_1_length_18919_cov_453.618_ID_122`))
pos2<-which(as.vector(fasta$AF086833.2)!=
as.vector(fasta$`C_Ebola-911_assemblypostmap_NODE_1_length_18919_cov_453.618_ID_122`))
pos2<-pos2[pos2>38 & pos2<18910]
# All positions variable either in 800 or in 911
pos<-sort(unique(c(pos1,pos2)))
# Print those
ref<-as.vector(fasta$AF086833.2)[pos]
names(ref)<-pos
ebola800<-as.vector(fasta$`Ebola-800_assemblypostmap_NODE_1_length_18872_cov_119.776_ID_37`)[pos]
ebola911<-as.vector(fasta$`C_Ebola-911_assemblypostmap_NODE_1_length_18919_cov_453.618_ID_122`)[pos]
ebola800<-as.vector(
fasta$`Ebola-800_assemblypostmap_NODE_1_length_18872_cov_119.776_ID_37`)[pos]
ebola911<-as.vector(
fasta$`C_Ebola-911_assemblypostmap_NODE_1_length_18919_cov_453.618_ID_122`)[pos]
var<-rbind(ref, ebola800, ebola911)
@
\item Print Variable positions between ref, Ebola-800 and Ebola-911:
\end{itemize}
<<>>=
var
@
Remarque: The alignment step involved an insertion pos 6918, so all the position > 6918 correspond to a position +1 in the alignemnt compared to the reference sequence alone.
Remarque: The alignment step involved an insertion pos 6918, so all the position $>$ 6918 correspond to a position +1 in the alignement compared to the reference sequence alone.
\section{Comparison of SNP calling with reference genome}
I will now analyse the variant identified using the reads directly aligned against the genome. To see if we recover the same positions?
Variants identified using the reads directly aligned against the genome will be analysed, to see if we recover the same positions.
\subsection{bcftool}
Thee first program used for this variant calling is bcftool.
\begin{itemize}
\item Read vcf (Variant Calling Files)
<<>>=
files<-paste0("../out_variant_trimmo/",
list.files("../out_variant_trimmo/", pattern="dedupbcftools.vcf$"))[1:2]
......@@ -128,9 +154,12 @@ parsevcf<-function(file){
print(".............................")
return(as.data.frame(locinfo[,c("POS")]))
}
@
\item Print variants:
\end{itemize}
<<>>=
resbcf<-sapply(files, parsevcf)
@
Only one snp is detected in Ebola-800 at position 1852 (T instead of C)
......@@ -140,19 +169,18 @@ Only 4 in ebola-911 (positions 1852, 2040 and 10557 and 10784). These 4 position
Freebayes is a more sophisticated variant calling tool.
<<echo=FALSE>>=
files<-paste0("../out_variant_trimmo/", list.files("../out_variant_trimmo/", pattern="ebola_dedup_freebayes.vcf$"))[1:2]
files<-paste0("../out_variant_trimmo/",
list.files("../out_variant_trimmo/",
pattern="ebola_dedup_freebayes.vcf$"))[1:2]
resfree<-sapply(files, parsevcf)
@
\section{Global summary}
A lot more de variant called. But they might not all be real? False positives?
<<fig.height=5>>=
\section{Global summary}
<<fig.height=5, fig.width=12>>=
plot(1:18959, rep(0, 18959), cex=0.1, xlab="position", ylab="", bty="n")
# graduation
points(seq(from=0, to=18000, by=1000), rep(0, 19), pch=3)
......@@ -171,14 +199,16 @@ points(tmp, rep(-0.2, length(tmp)), col="blue", cex=0.5)
tmp<-resfree$`../out_variant_trimmo/Ebola-911_S1_ebola_dedup_freebayes.vcf.locinfo[, c("POS")]`
points(tmp, rep(-0.3, length(tmp)), col="darkgreen", cex=0.5)
legend("topleft", c("genomes", "bcftool", "freebayes"), cex=1, col=c("red", "blue", "darkgreen"), pch=1, bty="n")
legend("topleft", c("genomes", "bcftool", "freebayes"),
cex=1,
col=c("red", "blue", "darkgreen"),
pch=1,
bty="n")
text(7000, 0.75, "Ebola-800", cex=1.5)
text(7000, -0.75, "Ebola-911", cex=1.5)
@
\end{document}
This is pdfTeX, Version 3.14159265-2.6-1.40.21 (TeX Live 2020) (preloaded format=pdflatex 2021.3.3) 4 MAR 2021 19:35
This is pdfTeX, Version 3.14159265-2.6-1.40.21 (TeX Live 2020) (preloaded format=pdflatex 2021.3.3) 9 MAR 2021 09:29
entering extended mode
restricted \write18 enabled.
%&-line parsing enabled.
......@@ -155,221 +155,222 @@ LaTeX Font Info: External font `cmex10' loaded for size
{/home/adminmarie/.TinyTeX/texmf-var/fonts/map/pdftex/updmap/pdftex.map}]
Package xcolor Warning: Incompatible color definition on input line 118.
Package xcolor Warning: Incompatible color definition on input line 127.
LaTeX Font Info: Font shape `OT1/cmtt/bx/n' in size <10.95> not available
(Font) Font shape `OT1/cmtt/m/n' tried instead on input line 120.
(Font) Font shape `OT1/cmtt/m/n' tried instead on input line 129.
Package xcolor Warning: Incompatible color definition on input line 132.
Package xcolor Warning: Incompatible color definition on input line 141.
Overfull \hbox (243.61221pt too wide) in paragraph at lines 134--134
[2]
Overfull \hbox (87.02927pt too wide) in paragraph at lines 143--143
[][]\OT1/cmtt/m/n/10.95 workdir[][]<-[][]"~/Documents/CIRI_BIBS_projects/2020_0
9_Reynard/seq_manual_aln/Ebola_aln_assemblyorigin3seq.fst"[][]
9_Reynard/seq_manual_aln/"[][]
[]
Overfull \hbox (352.83728pt too wide) in paragraph at lines 139--139
[][]\OT1/cmtt/m/n/10.95 pos1[][]<-[][]which[][]([][]as.vector[][](fasta[][]$[][
]AF086833.2)[][]!=[][]as.vector[][](fasta[][]$[][]^^REbola-800_assemblypostmap_
NODE_1_length_18872_cov_119.776_ID_37^^R))[][]
[]
Package xcolor Warning: Incompatible color definition on input line 151.
Overfull \hbox (234.01588pt too wide) in paragraph at lines 140--140
[][]\OT1/cmtt/m/it/10.95 #a few base missing at the start and at the end, they
are not significant cut between pos 38 and 18910[][]
Overfull \hbox (144.51614pt too wide) in paragraph at lines 155--155
[][]\OT1/cmtt/m/n/10.95 as.vector[][](fasta[][]$[][]^^REbola-800_assemblypostma
p_NODE_1_length_18872_cov_119.776_ID_37^^R))[][]
[]
Overfull \hbox (370.08334pt too wide) in paragraph at lines 144--144
[][]\OT1/cmtt/m/n/10.95 pos2[][]<-[][]which[][]([][]as.vector[][](fasta[][]$[][
]AF086833.2)[][]!=[][]as.vector[][](fasta[][]$[][]^^RC_Ebola-911_assemblypostma
p_NODE_1_length_18919_cov_453.618_ID_122^^R))[][]
Overfull \hbox (161.7622pt too wide) in paragraph at lines 162--162
[][]\OT1/cmtt/m/n/10.95 as.vector[][](fasta[][]$[][]^^RC_Ebola-911_assemblypost
map_NODE_1_length_18919_cov_453.618_ID_122^^R))[][]
[]
Overfull \hbox (197.62271pt too wide) in paragraph at lines 155--155
[][]\OT1/cmtt/m/n/10.95 ebola800[][]<-[][]as.vector[][](fasta[][]$[][]^^REbola-
800_assemblypostmap_NODE_1_length_18872_cov_119.776_ID_37^^R)[pos][][]
Overfull \hbox (133.01877pt too wide) in paragraph at lines 172--172
[] []\OT1/cmtt/m/n/10.95 fasta[][]$[][]^^REbola-800_assemblypostmap_NODE_1_l
ength_18872_cov_119.776_ID_37^^R)[pos][][]
[]
Overfull \hbox (214.86877pt too wide) in paragraph at lines 156--156
[][]\OT1/cmtt/m/n/10.95 ebola911[][]<-[][]as.vector[][](fasta[][]$[][]^^RC_Ebol
a-911_assemblypostmap_NODE_1_length_18919_cov_453.618_ID_122^^R)[pos][][]
Overfull \hbox (150.26483pt too wide) in paragraph at lines 174--174
[] []\OT1/cmtt/m/n/10.95 fasta[][]$[][]^^RC_Ebola-911_assemblypostmap_NODE_1
_length_18919_cov_453.618_ID_122^^R)[pos][][]
[]
Overfull \hbox (318.34515pt too wide) in paragraph at lines 167--167
Package xcolor Warning: Incompatible color definition on input line 184.
Overfull \hbox (53.90552pt too wide) in paragraph at lines 197--197
[]\OT1/cmtt/m/n/10.95 ## 44 45 48 928 1852 2040 2410 2411 2551 2762
6182 6185 6918 6949 9555 10558 10785 10905 14188 18139 18454[]
6182 6185 6918[]
[]
Overfull \hbox (306.84778pt too wide) in paragraph at lines 167--167
Overfull \hbox (48.15683pt too wide) in paragraph at lines 197--197
[]\OT1/cmtt/m/n/10.95 ## ref "g" "c" "a" "t" "c" "c" "t" "t" "g" "a"
"c" "c" "-" "a" "c" "g" "t" "c" "g" "g" "c"[]
"c" "c" "-"[]
[]
Overfull \hbox (306.84778pt too wide) in paragraph at lines 167--167
Overfull \hbox (48.15683pt too wide) in paragraph at lines 197--197
[]\OT1/cmtt/m/n/10.95 ## ebola800 "t" "t" "g" "c" "t" "t" "t" "t" "a" "c"
"a" "t" "a" "a" "t" "a" "c" "t" "a" "a" "t"[]
"a" "t" "a"[]
[]
Overfull \hbox (306.84778pt too wide) in paragraph at lines 167--167
Overfull \hbox (48.15683pt too wide) in paragraph at lines 197--197
[]\OT1/cmtt/m/n/10.95 ## ebola911 "g" "c" "a" "c" "t" "t" "c" "c" "a" "c"
"a" "t" "-" "c" "t" "a" "c" "t" "a" "a" "c"[]
"a" "t" "-"[]
[]
[2]
Overfull \hbox (11.07118pt too wide) in paragraph at lines 172--172
[3]
Overfull \hbox (11.07118pt too wide) in paragraph at lines 203--203
[]\OT1/cmr/bx/n/14.4 Comparison of SNP call-ing with ref-er-ence genome
[]
Package xcolor Warning: Incompatible color definition on input line 179.
Package xcolor Warning: Incompatible color definition on input line 215.
[3]
Overfull \hbox (134.38715pt too wide) in paragraph at lines 182--182
Overfull \hbox (161.7622pt too wide) in paragraph at lines 218--218
[] []\OT1/cmtt/m/n/10.95 list.files[][]([][]"../out_variant_trimmo
/"[][],[] []pattern[][]=[][]"dedupbcftools.vcf$"[][]))[[][]1[][]:[][]2[][]][][]
[]
LaTeX Font Info: Trying to load font information for OMS+cmtt on input line
184.
220.
(/usr/share/R/share/texmf/tex/latex/omscmtt.fd
File: omscmtt.fd
)
LaTeX Font Info: Font shape `OMS/cmtt/m/n' in size <10.95> not available
(Font) Font shape `OMS/cmsy/m/n' tried instead on input line 184.
(Font) Font shape `OMS/cmsy/m/n' tried instead on input line 220.
Overfull \hbox (25.16208pt too wide) in paragraph at lines 193--193
Overfull \hbox (0.79895pt too wide) in paragraph at lines 227--227
[] []\OT1/cmtt/m/n/10.95 locinfo[][]<-[][]locinfo[[][]selectQUAL[][](locinfo
,[] []threshold[] []=[] []20[][]),][][]
[]
Overfull \hbox (52.53714pt too wide) in paragraph at lines 229--229
[] []\OT1/cmtt/m/n/10.95 print[][]([][]as.data.frame[][](locinfo[,[][]c[][](
[][]"CHROM"[][],[] []"POS"[][],[] []"REF"[][],[] []"ALT"[][])]))[][]
[]
Overfull \hbox (71.15158pt too wide) in paragraph at lines 222--222
Package xcolor Warning: Incompatible color definition on input line 240.
[4]
Overfull \hbox (71.15158pt too wide) in paragraph at lines 266--266
[]\OT1/cmtt/m/n/10.95 ## Scanning file ../out_variant_trimmo/Ebola-800_S2_ebola
_dedupbcftools.vcf[]
[]
Overfull \hbox (71.15158pt too wide) in paragraph at lines 222--222
Overfull \hbox (71.15158pt too wide) in paragraph at lines 266--266
[]\OT1/cmtt/m/n/10.95 ## Scanning file ../out_variant_trimmo/Ebola-911_S1_ebola
_dedupbcftools.vcf[]
[]
[4]
Package xcolor Warning: Incompatible color definition on input line 235.
Package xcolor Warning: Incompatible color definition on input line 278.
Overfull \hbox (82.64896pt too wide) in paragraph at lines 312--312
Overfull \hbox (82.64896pt too wide) in paragraph at lines 355--355
[]\OT1/cmtt/m/n/10.95 ## Scanning file ../out_variant_trimmo/Ebola-800_S2_ebola
_dedup_freebayes.vcf[]
[]
Overfull \hbox (82.64896pt too wide) in paragraph at lines 312--312
Overfull \hbox (82.64896pt too wide) in paragraph at lines 355--355
[]\OT1/cmtt/m/n/10.95 ## Scanning file ../out_variant_trimmo/Ebola-911_S1_ebola
_dedup_freebayes.vcf[]
[]
[5]
Underfull \vbox (badness 10000) detected at line 313
Underfull \vbox (badness 10000) detected at line 356
[]
Underfull \vbox (badness 10000) detected at line 313
Underfull \vbox (badness 10000) detected at line 356
[]
[6]
Package xcolor Warning: Incompatible color definition on input line 320.
Package xcolor Warning: Incompatible color definition on input line 364.
Overfull \hbox (53.90552pt too wide) in paragraph at lines 322--322
[7]
Overfull \hbox (53.90552pt too wide) in paragraph at lines 366--366
[][]\OT1/cmtt/m/n/10.95 plot[][]([][]1[][]:[][]18959[][],[] []rep[][]([][]0[][]
,[] []18959[][]),[] []cex[][]=[][]0.1[][],[] []xlab[][]=[][]"position"[][],[] [
]ylab[][]=[][]""[][],[] []bty[][]=[][]"n"[][])[][]
[]
Overfull \hbox (2.16733pt too wide) in paragraph at lines 325--325
Overfull \hbox (2.16733pt too wide) in paragraph at lines 369--369
[][]\OT1/cmtt/m/n/10.95 points[][]([][]seq[][]([][]from[][]=[][]0[][],[] []to[]
[]=[][]15000[][],[] []by[][]=[][]5000[][]),[] []rep[][]([][]0[][],[] []4[][]),[
] []pch[][]=[][]3[][],[] []cex[][]=[][]3[][])[][]
[]
Overfull \hbox (168.87927pt too wide) in paragraph at lines 328--328
Overfull \hbox (168.87927pt too wide) in paragraph at lines 372--372
[][]\OT1/cmtt/m/n/10.95 tmp[][]<-[][]resbcf[][]$[][]^^R../out_variant_trimmo/Eb
ola-800_S2_ebola_dedupbcftools.vcf.locinfo[, c("POS")]^^R[][]
[]
Overfull \hbox (186.12534pt too wide) in paragraph at lines 330--330
Overfull \hbox (186.12534pt too wide) in paragraph at lines 374--374
[][]\OT1/cmtt/m/n/10.95 tmp[][]<-[][]resfree[][]$[][]^^R../out_variant_trimmo/E
bola-800_S2_ebola_dedup_freebayes.vcf.locinfo[, c("POS")]^^R[][]
[]
Overfull \hbox (168.87927pt too wide) in paragraph at lines 335--335
Overfull \hbox (168.87927pt too wide) in paragraph at lines 379--379
[][]\OT1/cmtt/m/n/10.95 tmp[][]<-[][]resbcf[][]$[][]^^R../out_variant_trimmo/Eb
ola-911_S1_ebola_dedupbcftools.vcf.locinfo[, c("POS")]^^R[][]
[]
Overfull \hbox (186.12534pt too wide) in paragraph at lines 337--337
Overfull \hbox (186.12534pt too wide) in paragraph at lines 381--381
[][]\OT1/cmtt/m/n/10.95 tmp[][]<-[][]resfree[][]$[][]^^R../out_variant_trimmo/E
bola-911_S1_ebola_dedup_freebayes.vcf.locinfo[, c("POS")]^^R[][]
[]
<figure/unnamed-chunk-7-1.pdf, id=36, 867.24pt x 361.35pt>
File: figure/unnamed-chunk-7-1.pdf Graphic file (type pdf)
<use figure/unnamed-chunk-7-1.pdf>
Package pdftex.def Info: figure/unnamed-chunk-7-1.pdf used on input line 394.
(pdftex.def) Requested size: 360.0pt x 150.00156pt.
Overfull \hbox (289.60172pt too wide) in paragraph at lines 340--340
[][]\OT1/cmtt/m/n/10.95 legend[][]([][]"topleft"[][],[] []c[][]([][]"genomes"[]
[],[] []"bcftool"[][],[] []"freebayes"[][]),[] []cex[][]=[][]1[][],[] []col[][]
=[][]c[][]([][]"red"[][],[] []"blue"[][],[] []"darkgreen"[][]),[] []pch[][]=[][
]1[][],[] []bty[][]=[][]"n"[][])[][]
[]
<figure/unnamed-chunk-4-1.pdf, id=29, 505.89pt x 361.35pt>
File: figure/unnamed-chunk-4-1.pdf Graphic file (type pdf)
<use figure/unnamed-chunk-4-1.pdf>
Package pdftex.def Info: figure/unnamed-chunk-4-1.pdf used on input line 346.
(pdftex.def) Requested size: 360.0pt x 257.14474pt.
Overfull \hbox (17.0pt too wide) in paragraph at lines 346--347
Overfull \hbox (17.0pt too wide) in paragraph at lines 394--395
[][]
[]
[7] [8 <./figure/unnamed-chunk-4-1.pdf>] (./2_ebola_report.aux) )
[8 <./figure/unnamed-chunk-7-1.pdf>] (./2_ebola_report.aux) )
Here is how much of TeX's memory you used:
2982 strings out of 480971
43357 string characters out of 5907380
356336 words of memory out of 5000000
20188 multiletter control sequences out of 15000+600000
408937 words of font info for 46 fonts, out of 8000000 for 9000
2987 strings out of 480971
43420 string characters out of 5907380
356599 words of memory out of 5000000
20191 multiletter control sequences out of 15000+600000
409624 words of font info for 48 fonts, out of 8000000 for 9000
14 hyphenation exceptions out of 8191
68i,6n,74p,434b,327s stack positions out of 5000i,500n,10000p,200000b,80000s
</home/adminma
rie/.TinyTeX/texmf-dist/fonts/type1/public/amsfonts/cm/cmbx10.pfb></home/adminm
arie/.TinyTeX/texmf-dist/fonts/type1/public/amsfonts/cm/cmbx12.pfb></home/admin
marie/.TinyTeX/texmf-dist/fonts/type1/public/amsfonts/cm/cmitt10.pfb></home/adm
inmarie/.TinyTeX/texmf-dist/fonts/type1/public/amsfonts/cm/cmr10.pfb></home/adm
inmarie/.TinyTeX/texmf-dist/fonts/type1/public/amsfonts/cm/cmr12.pfb></home/adm
inmarie/.TinyTeX/texmf-dist/fonts/type1/public/amsfonts/cm/cmr17.pfb></home/adm
inmarie/.TinyTeX/texmf-dist/fonts/type1/public/amsfonts/cm/cmsy10.pfb></home/ad
68i,6n,74p,352b,388s stack positions out of 5000i,500n,10000p,200000b,80000s
</home/adminmarie
/.TinyTeX/texmf-var/fonts/pk/ljfour/jknappen/ec/tcrm1095.600pk></home/adminmari
e/.TinyTeX/texmf-dist/fonts/type1/public/amsfonts/cm/cmbx10.pfb></home/adminmar
ie/.TinyTeX/texmf-dist/fonts/type1/public/amsfonts/cm/cmbx12.pfb></home/adminma
rie/.TinyTeX/texmf-dist/fonts/type1/public/amsfonts/cm/cmitt10.pfb></home/admin
marie/.TinyTeX/texmf-dist/fonts/type1/public/amsfonts/cm/cmmi10.pfb></home/admi
nmarie/.TinyTeX/texmf-dist/fonts/type1/public/amsfonts/cm/cmr10.pfb></home/admi
nmarie/.TinyTeX/texmf-dist/fonts/type1/public/amsfonts/cm/cmr12.pfb></home/admi
nmarie/.TinyTeX/texmf-dist/fonts/type1/public/amsfonts/cm/cmr17.pfb></home/admi
nmarie/.TinyTeX/texmf-dist/fonts/type1/public/amsfonts/cm/cmsy10.pfb></home/adm
inmarie/.TinyTeX/texmf-dist/fonts/type1/public/amsfonts/cm/cmti10.pfb></home/ad
minmarie/.TinyTeX/texmf-dist/fonts/type1/public/amsfonts/cm/cmtt10.pfb>
Output written on 2_ebola_report.pdf (8 pages, 174971 bytes).
Output written on 2_ebola_report.pdf (8 pages, 190395 bytes).
PDF statistics:
70 PDF objects out of 1000 (max. 8388607)
49 compressed objects within 1 object stream
83 PDF objects out of 1000 (max. 8388607)
59 compressed objects within 1 object stream
0 named destinations out of 1000 (max. 500000)
6 words of extra memory for PDF output out of 10000 (max. 10000000)
No preview for this file type
......@@ -98,22 +98,31 @@
\section{Introduction}
For detailed summary of the analysis, see README.md file.
For detailed summary of the analysis, see \textit{README.md} file.
~
Briefly, illumina paired-end reads from 2 Ebola strains, Ebola-800 and Ebola-911 were analysed.
The aim of this study is to identify all differences between the sequences of 800 and 911 compared to the reference sequences.
For this, 2 approaches are possible, first, to assemble genomes de novo and the compare genome sequences to reference genome sequence, second, to align reads against the reference genomes and use software which identify variant directly from the read alignement.
After quality check and read cleaning by Trimmomatic, I 1/ assembled both genome de novo using SPADES 2/ mapped read against a reference genome (Bowtie) and called variants using 2 softwares (bcftools and Freebayes) 3/ Reassembled reads using only mapped reads (SPADES).
~
The second assembly using only reads that mapped against the reference genome was of higher quality. I will use this one in this analysis.
~
The aim of this report is to identify and compare the variations found using these different approaches.
\section{Comparison of genome assemblies with reference genome}
Sequences from the reference genome and the 2 SPADES assembly (2nd version, mapped reads) were manually aligned. See report on the GP gene Ebola\_assembly.pdf.
~
\begin{itemize}
\item Load R packages
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
......@@ -127,54 +136,81 @@ Sequences from the reference genome and the 2 SPADES assembly (2nd version, mapp
\end{kframe}
\end{knitrout}
Read the alignment.
\item Read the alignment.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{workdir}\hlkwb{<-}\hlstr{"~/Documents/CIRI_BIBS_projects/2020_09_Reynard/seq_manual_aln/"}
\hlstd{file}\hlkwb{<-}\hlkwd{paste0}\hlstd{(workdir,} \hlstr{"Ebola_aln_assemblyorigin3seq.fst"}\hlstd{)}
\hlstd{fasta}\hlkwb{<-}\hlkwd{read.fasta}\hlstd{(file)}
\end{alltt}
\end{kframe}
\end{knitrout}
\item Pairwise comparison: Which are the variable positions
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{workdir}\hlkwb{<-}\hlstr{"~/Documents/CIRI_BIBS_projects/2020_09_Reynard/seq_manual_aln/Ebola_aln_assemblyorigin3seq.fst"}
\hlcom{#fasta<-read.dna(workdir, format = "fasta", as.matrix = TRUE)}
\hlstd{fasta}\hlkwb{<-}\hlkwd{read.fasta}\hlstd{(workdir)}
\hlcom{# variable positions between ref and 800}
\hlstd{pos1}\hlkwb{<-}\hlkwd{which}\hlstd{(}\hlkwd{as.vector}\hlstd{(fasta}\hlopt{$}\hlstd{AF086833.2)}\hlopt{!=}\hlkwd{as.vector}\hlstd{(fasta}\hlopt{$}\hlstd{`Ebola-800_assemblypostmap_NODE_1_length_18872_cov_119.776_ID_37`))}
\hlcom{#a few base missing at the start and at the end, they are not significant cut between pos 38 and 18910}
\hlstd{pos1}\hlkwb{<-}\hlkwd{which}\hlstd{(}\hlkwd{as.vector}\hlstd{(fasta}\hlopt{$}\hlstd{AF086833.2)}\hlopt{!=}
\hlkwd{as.vector}\hlstd{(fasta}\hlopt{$}\hlstd{`Ebola-800_assemblypostmap_NODE_1_length_18872_cov_119.776_ID_37`))}
\hlcom{# a few base missing at the start and at the end, }
\hlcom{# they are not significant cut between pos 38 and 18910}
\hlstd{pos1}\hlkwb{<-}\hlstd{pos1[pos1}\hlopt{>}\hlnum{38} \hlopt{&} \hlstd{pos1}\hlopt{<}\hlnum{18910}\hlstd{]}
\hlcom{# variable positions between ref and 911}
\hlstd{pos2}\hlkwb{<-}\hlkwd{which}\hlstd{(}\hlkwd{as.vector}\hlstd{(fasta}\hlopt{$}\hlstd{AF086833.2)}\hlopt{!=}\hlkwd{as.vector}\hlstd{(fasta}\hlopt{$}\hlstd{`C_Ebola-911_assemblypostmap_NODE_1_length_18919_cov_453.618_ID_122`))}
\hlstd{pos2}\hlkwb{<-}\hlkwd{which}\hlstd{(}\hlkwd{as.vector}\hlstd{(fasta}\hlopt{$}\hlstd{AF086833.2)}\hlopt{!=}
\hlkwd{as.vector}\hlstd{(fasta}\hlopt{$}\hlstd{`C_Ebola-911_assemblypostmap_NODE_1_length_18919_cov_453.618_ID_122`))}
\hlstd{pos2}\hlkwb{<-}\hlstd{pos2[pos2}\hlopt{>}\hlnum{38} \hlopt{&} \hlstd{pos2}\hlopt{<}\hlnum{18910}\hlstd{]}
\hlcom{# All positions variable either in 800 or in 911}
\hlstd{pos}\hlkwb{<-}\hlkwd{sort}\hlstd{(}\hlkwd{unique}\hlstd{(}\hlkwd{c}\hlstd{(pos1,pos2)))}
\hlcom{# Print those}
\hlstd{ref}\hlkwb{<-}\hlkwd{as.vector}\hlstd{(fasta}\hlopt{$}\hlstd{AF086833.2)[pos]}
\hlkwd{names}\hlstd{(ref)}\hlkwb{<-}\hlstd{pos}
\hlstd{ebola800}\hlkwb{<-}\hlkwd{as.vector}\hlstd{(fasta}\hlopt{$}\hlstd{`Ebola-800_assemblypostmap_NODE_1_length_18872_cov_119.776_ID_37`)[pos]}
\hlstd{ebola911}\hlkwb{<-}\hlkwd{as.vector}\hlstd{(fasta}\hlopt{$}\hlstd{`C_Ebola-911_assemblypostmap_NODE_1_length_18919_cov_453.618_ID_122`)[pos]}
\hlstd{ebola800}\hlkwb{<-}\hlkwd{as.vector}\hlstd{(}
\hlstd{fasta}\hlopt{$}\hlstd{`Ebola-800_assemblypostmap_NODE_1_length_18872_cov_119.776_ID_37`)[pos]}
\hlstd{ebola911}\hlkwb{<-}\hlkwd{as.vector}\hlstd{(}
\hlstd{fasta}\hlopt{$}\hlstd{`C_Ebola-911_assemblypostmap_NODE_1_length_18919_cov_453.618_ID_122`)[pos]}
\hlstd{var}\hlkwb{<-}\hlkwd{rbind}\hlstd{(ref, ebola800, ebola911)}
\end{alltt}
\end{kframe}
\end{knitrout}
\item Print Variable positions between ref, Ebola-800 and Ebola-911:
\end{itemize}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{var}
\end{alltt}
\begin{verbatim}
## 44 45 48 928 1852 2040 2410 2411 2551 2762 6182 6185 6918 6949 9555 10558 10785 10905 14188 18139 18454
## ref "g" "c" "a" "t" "c" "c" "t" "t" "g" "a" "c" "c" "-" "a" "c" "g" "t" "c" "g" "g" "c"
## ebola800 "t" "t" "g" "c" "t" "t" "t" "t" "a" "c" "a" "t" "a" "a" "t" "a" "c" "t" "a" "a" "t"
## ebola911 "g" "c" "a" "c" "t" "t" "c" "c" "a" "c" "a" "t" "-" "c" "t" "a" "c" "t" "a" "a" "c"
## 44 45 48 928 1852 2040 2410 2411 2551 2762 6182 6185 6918
## ref "g" "c" "a" "t" "c" "c" "t" "t" "g" "a" "c" "c" "-"
## ebola800 "t" "t" "g" "c" "t" "t" "t" "t" "a" "c" "a" "t" "a"
## ebola911 "g" "c" "a" "c" "t" "t" "c" "c" "a" "c" "a" "t" "-"
## 6949 9555 10558 10785 10905 14188 18139 18454
## ref "a" "c" "g" "t" "c" "g" "g" "c"
## ebola800 "a" "t" "a" "c" "t" "a" "a" "t"
## ebola911 "c" "t" "a" "c" "t" "a" "a" "c"
\end{verbatim}
\end{kframe}
\end{knitrout}
Remarque: The alignment step involved an insertion pos 6918, so all the position > 6918 correspond to a position +1 in the alignemnt compared to the reference sequence alone.
Remarque: The alignment step involved an insertion pos 6918, so all the position $>$ 6918 correspond to a position +1 in the alignement compared to the reference sequence alone.
\section{Comparison of SNP calling with reference genome}
I will now analyse the variant identified using the reads directly aligned against the genome. To see if we recover the same positions?
Variants identified using the reads directly aligned against the genome will be analysed, to see if we recover the same positions.
\subsection{bcftool}
Thee first program used for this variant calling is bcftool.
\begin{itemize}
\item Read vcf (Variant Calling Files)
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
......@@ -194,7 +230,15 @@ I will now analyse the variant identified using the reads directly aligned again
\hlkwd{print}\hlstd{(}\hlstr{"............................."}\hlstd{)}
\hlkwd{return}\hlstd{(}\hlkwd{as.data.frame}\hlstd{(locinfo[,}\hlkwd{c}\hlstd{(}\hlstr{"POS"}\hlstd{)]))}
\hlstd{\}}
\end{alltt}
\end{kframe}
\end{knitrout}