原文地址:de novo assembly exercises(转载)作者:珂珂 de novo assembly exercises simon rasm ssen - j ne 2013 overview in this exercise we will try to perform de novo assembly of ill mina paired end reads. the data is from a vibrio cholerae strain isolated in nepal. yo will try to: r n fastqc, adaptor and q ality trimming reads co nt k-mers and estimate genome size correct reads sing q ake determine insert size of paired end reads r n de novo assembly sing soapdenovo calc late assembly statistics plot coverage and length histograms of the assembly re-assembly to improve assembly vis alize assembly sing circoletto fastqc and trimming create dir and copy data cd homelocalngs_co rsest dxxx mkdir denovo cd denovo cp homelocalexercisesdenovovchol* . start by r nning fastqc mkdir fastqc fastqc -o fastqc *.txt.gz firefox fastqcvchol-001_6_1_seq ence_fastqcfastqc_reportl & firefox fastqcvchol-001_6_2_seq ence_fastqcfastqc_reportl & there are q ite some iss es with this data (b t dont spend to m ch time on st dying the report), we sho ld clean it p first. lets find o t what type of encoding the q alities are in: gzcat vchol-001_6_1_seq ence.txt.gz | fastx_detect_fq.py q1. which format are they encoded in? lets trim the reads sing c tadapt. in the command below the most freq ent adapterprimer seq ence is already pasted in - dont worry abo t the others in o r case they are j st variations of that one. c tadapt will ass me a 10% error rate so they sho ld also be clipped as well. also we we se minim m of 40nt and trim to q ality 20 (~5 mins). note that both commands ends with "&" meaning that they will both r n on the same time in the backgro nd. yo m st wait for both commands to finish before going to the next command! c tadapt -b gatcggaagagcacacgtctgaactccagtcacatcacgatatcgtatgc -o 3 -m 40 -q 20 --q ality-base=64 vchol-001_6_1_seq ence.txt.gz  vchol-001_6_1.c t.fastq & c tadapt -b gatcggaagagcgtcgtgtagggaaagagggtagatctcggtggtcgccg -o 3 -m 40 -q 20 --q ality-base=64 vchol-001_6_2_seq ence.txt.gz  vchol-001_6_2.c t.fastq & then we need to p t the reads in the same order sing cmpfastq (1min)  vchol-001_6_1.c t.fastq vchol-001_6_2.c t.fastq ok, now we have o r reads trimmed and ready to go. co nt stats of the trimmed data s ch as avg readlength, min length, max length, no. reads and no. bases. note down the avg readlength and the no. of bases. fastx_readlength.py --i vchol-001_6_1.c t.fastq-common.o t fastx_readlength.py --i vchol-001_6_2.c t.fastq-common.o t genome size estimation lets try to co nt the occ rence of k-mers in the data - a k-mer is simply a string of n cleotides of a certain length. lets say we co nt all the 15-mers that are in the read data we have - we do this sing a program called jellyfish. here we tell jellyfish to co nt 15-mers and also add co nts from the complementary strand. afterwards we tell it to create a histogram that we will plot sing r jellyfish co nt -t 2 -m 15 -s  -o vchol-001 -c vchol-001_6_*.c t.fastq-common.o t jellyfish histo vchol-001_0  vchol-001.histo r dat=read.table("vchol-001.histo") barplot(dat[,2], xlim=c(0,150), ylim=c(0,5e5), ylab="no of kmers", xlab="co nts of a k-mer", names.arg=dat[,1], cex.names=0.8) dev.print("vchol-001.histo.pdf", device=pdf) the plot shows how many times a k-mer is fo nd (x-axis) and the n mber of kmers with this co nt (y-axis). the bars in the lower end are the k-mers that only occ r very few times, these are probably seq encing errors, whereas the k-mers that occ rs many times are "real" k-mers. we can se a program like q ake to se the information from the "real" k-mers to correct similar "error k-mers". in this way we can increase the performance of o r assembly (this also works for snp calling as well). q2. where is the peak at? before we correct the reads, we will try to estimate the genome size of o r genome. of co rse we know that it is a v. cholerae, b t if we did not know that or it was an nknown species then we co ld do like this. n = (m*l)(l-k+1) and genome_size = tn, where n: depth, m: kmer peak, k: kmer-size, l: avg readlength, t: total bases. try to p t in the n mbers (some of them yo have from "fastx_readlength.py" command) and see what yo get. the reference vibrio cholerae genome is aro nd 4mb, yo sho ld get within +- 10% of this. q3. what is the estimated genome size q ake error correction lets try to correct the reads sing q ake. it is act ally rather easy to do, we need to make a file with the read data that we have and then start q ake with k=15 and tell it that the q alities are ill mina (-q 64). the corrected reads will be *.cor.o t # this takes some time - you can copy my files # echo "vchol-001_6_1.c t.fastq-common.o t vchol-001_6_2.c t.fastq-common.o t"  file_list homelocalbinq akebinq  -f file_list -k 15 -p 2 -q 64 --header # copy my files # ln -s homelocalngs_co rsest d100denovo*.cor.o t . de novo assembly now the reads are ready for de novo assembly. we are going to se soapdenovo as the assembler, it ses the de br ijn approach. when r nning denovo assemblies one sho ld try different k-mer sizes - the k-mer size is what is being sed for b ilding the de br ijn graph and is therefore very important. there is c rrently no way to estimate what will be the optimal k before r nning, the best k-size may change between different dataset and may change between different assemblers for the same dataset. soapdenovo req ires a config ration file with different information, a sample file is the vchol-001.soap.conf. yo need to open it and change the path to yo r file. to open it write "n vchol-001.soap.conf" and then navigate to "q1" and "q2", here paste in the path and filename of yo r two files. note that we have written "avg_ins=200" even tho gh we do not know the average insert size. we will r n an initial assembly to create contigs and then extract the information and then rer n the f ll assembly again. when yo have changed the paths save the file (ctrl-s) then start the assembler. we will try first with a k of 35. soapdenovo-127mer pregraph -s vchol-001.soap.conf -k 35 -p 2 -o initial soapdenovo-127mer contig -g initial now yo sho ld have a file called "initial.contig", we need to map o r reads back to the contigs to identify the insert size, j st as we did in the alignment exercise. lets only map the first 100.000 reads - this sho ld be eno gh. head -n  vchol-001_6_1.c t.fastq-common.o t  vchol_sample_1.fastq head -n  vchol-001_6_2.c t.fastq-common.o t  vchol_sample_2.fastq bwa index initial.contig bwa mem initial.contig vchol_sample_1.fastq vchol_sample_2.fastq | samtools view -sb -  initial.sample.bam samtools view initial.sample.bam | c t -f9  initial.insertsizes.txt r a = read.table("initial.insertsizes.txt") a.v = a[a[,1]0,1] mn = q antile(a.v, seq(0,1,0.05))[4] mx = q antile(a.v, seq(0,1,0.05))[18] mean(a.v[a.v = mn & a.v = mx]) # mean sd(a.v[a.v = mn & a.v = mx]) # sd q4. what are the mean insert size and standard deviation of o r library? open "vchol-001.soap.conf" and change the avg_ins to the val e yo fo nd. then we are going to rer n the assembly where we also do scaffolding. beca se the assembly takes some time to finish, yo will all chose a different k-mer to r n and then we will compare to find the "best" assembly. chose one of the k-sizes on the google doc ment and change the "xx" val es below to yo r n mber. remember to write yo r name and save the doc ment. soapdenovo-127mer all -s vchol-001.soap.conf -k xx -p 2 -o asmkxx when it is done, we need to calc late some statistics on it, there is a script from the assemblathon that does it nicely, b t first lets remove all seq ences shorter than 100 fastx_filterfasta.py --i asmkxx.scafseq --min 100 assemblathon_stats.pl asmkxx.scafseq.filtered_100.fa  asmkxx.scafseq.stats open the asmkxx.scafseq.stats file and fill in the information in the google doc ment for yo r respective k. we will then compare the res lts and find the best k-size for the partic lar dataset. if there are still k-sizes that has not been r n, feel free to r n them. copy the best assembly to yo r folder and se this assembly from now on (we renamed it to vchol-001.best.fa). cp homelocalexercisesdenovobestvchol-001.best.fa . cp homelocalexercisesdenovobestvchol-001.best.stats . q5. how does the s m of the assembly compare with what we estimated sing kmers? q6. how many contigs were joined together in scaffolds? why do yo think it is a s ch a low n mber (hint: repeats, pe insert size)? coverage of assembly lets calc late coverage and length for each seq ence and plot it as a histogram in r. we also plot the lengths of the scaffolds. fastx_soapcov.py --i vchol-001.best.fa  vchol-001.best.cov r library(plotrix) dat=read.table("vchol-001.best.cov", sep="t") par(mfrow=c(1,2)) weighted.hist(w=dat[,2], x=dat[,1], breaks=seq(0,100, 1), main="weighted coverage", xlab="contig coverage") hist(dat[,1], xlim=c(0,100), breaks=seq(0,1000,1), main="raw coverage", xlab="contig coverage") dev.print("best.coverage.pdf", device=pdf) # lengths par(mfrow=c(1,1)) barplot(rev(sort(dat[,2])), xlab="# scaffold", ylab="length", main="scaffold lengths") dev.print("scaffold.lengths.pdf", device=pdf) acroread best.coverage.pdf & acroread scaffold.lengths.pdf & the left plot shows the length-weighted coverage (of k-mers) of the contigsscaffolds - this means that long seq ences has a larger weight compared to short seq ences. the plot on the right side shows a normal histogram of contig coverage. by comparing the two plots we see that the majority of the assembly has a k-mer coverage aro nd 20-30x, and that remaining assembly are short seq ences. looking at the plot with the lengths yo see that the majority of the assembly are in q ite long scaffolds. q7. can yo think of why some short contigs have m ch higher coverage compared to main assembly? q8. can yo think of why some short contigs have m ch lower coverage compared to the main assembly? re-assembly now we created a fairly good assembly, b t lets see if we can do it better. lets try to map the reads to the assembly and then only se mapped reads for another assembly. we map sing bwa. bwa index vchol-001.best.fa bwa mem -m -t 4 vchol-001.best.fa vchol-001_6_1.c t.fastq-common.cor.o t vchol-001_6_2.c t.fastq-common.cor.o t | samtools view -sb -  vchol-001.bwa.bam when it is done we r n we look at the statistics of the assembly. samtools flagstat vchol-001.bwa.bam q9. how many reads are properly paired (properly paired means that the map within expected range of each other)? lets extract the properly paired reads from the alignment and make a new assembly sing only these. removing reads from the assembly that may make tro bles will decrease the complexity of the de br ijn graph and hopef lly lead to a better assembly with less errors. go to explain-flagsand find the flag-n mber for a properly paired read. insert the n mber below instead of y to create a bam-file of only properly paired reads and the extract to fastq-files. q10. what is the flag for properly paired? samtools view -b -f y vchol-001.bwa.bam  vchol-001.bwa.pp.bam java -jar homelocalbinsamtofastq.jar input=vchol-001.bwa.pp.bam fastq=vchol-001_1.remade.fastq second_end_fastq=vchol-001_2.remade.fastq update the "vchol-001.soap.conf" file with the filenames "vchol-001_1.remade.fastq" and "vchol-001_2.remade.fastq" and then r n soapdenovo again. now we j st se k=73 only, filter the o tp t scaffolds and r n stats on it. soapdenovo-127mer all -s vchol-001.soap.conf -k 73 -p 2 -o vchol-001.reasm fastx_filterfasta.py --i vchol-001.reasm.scafseq --min 100 assemblathon_stats.pl vchol-001.reasm.scafseq.filtered_100.fa  vchol-001.reasm.scafseq.filtered_100.stats open the new "vchol-001.reasm.scafseq.filtered_100.stats" and the original file "vchol-001.best.stats" and compare the res lts. q11. how m ch improvement is there in the assembly? is this what yo expect from the n mber of properly paired reads in q9? now that we are done with the assemblies lets clean p and remove o r reads (we are sing 7gb of disk-space) rm *.fastq rm *.o t vis alization sing circoletto lets try to vis alize the assembly. beca se we know that the data is from a vibrio cholerae we can try to compare o r assembly with the a v.cholerae reference genome. first lets filter the assembly to minim m 500bp. fastx_filterfasta.py --i vchol-001.reasm.scafseq --min 500 open a browser on yo r own comp ter and go to the circoletto webpage. circoletto is an easy to se web-program that b ilds pon circos - a program to create amazing plots with genomic (and other) data. open the filtered assembly ("n vchol-001.reasm.scafseq.filtered_500.fa") select all and copy-paste it into the pper box (q ery fasta). do the same with the reference genome ("n vibrio_cholerae_o1_n.fa") and copy-paste it into the box j st below. in the "o tp t" section click the "only show the best hit per q ery" and click "s bmit to circoletto". when it has finished processing click "see the circoletto o tp t". if it is not working yo can see the file here. yo sho ld see the two chromosomes of v. cholerae on the left hand side of the fig re (named "gi|...") and then the alignment of o r assembly to it. the color of the alignments are the bitscore, the red are the most confident and the black the least confident. q12. do yo think that o r genome is similar to the reference genome? q13. are there scaffoldscontigs that does not or only partially map to the reference genome? q14. what do yo think the region on chromosome 2 (the small one) that has many small low-confident mappings co ld be. hint: see the v. cholerae genome paper and search for "v. cholerae integron island"? bon s assembly st ff if yo are f rther interested in checking yo r assembly, it is possible to call snps by remapping reads to the assembly. for haploid chromosomes we do not expect any snps, so therefore any snps that are called are markers of assembly error or problems. additionally one may also se re-mapping of paired end reads to break contigs where there are no s pport from read pairs. assembly of 454 or ion torrent data (optional) for 454 and ion torrent data we can se newbler, an assembler made by roche the company prod cing the 454 machines. the fact that the error-profile of ion torrent data is similar to 454 data, we can get q ite good assemblies sing newbler for ion data. the commands to r n newbler assembly is this. it can take inp t as fasta+q al, fasta, sff and fastq. note, that it if ion torrent data is in fastq-format it m st be converted to fasta, beca se else newbler will think it is ill mina reads and prod ce a very bad assembly. setenv path homelocalbin454bin:${path} newassembly o tdir addr n -lib libname o tdir single.fasta addr n -p -lib pe1 o tdir paired_1.fasta addr n -p -lib pe1 o tdir paired_2.fasta r nproject -cp 2 o tdir lets try with the e.coli data from the german c c mberspro t o tbreak. mkdir ecoli cd ecoli cp homelocalexercisesdenovoecoli* . 454 data as fastq: setenv path homelocalbin454bin:${path} newassembly asm454 addr n -lib wgs asm454 gos1.fastq r nproject -cp 2 asm454 ion torrent data. first we need to convert the fastq to fasta. setenv path homelocalbin454bin:${path} fastx_fastq2fasta.py --i iontorrent.fq --o iontorrent newassembly asmion addr n -lib wgs asmion iontorrent.fasta r nproject -cp 2 asmion calc late statstics of the assemblies and try to compare them. assemblathon_stats.pl asm454assembly454largecontigs.fna  asm454.largecontigs.stats assemblathon_stats.pl asmionassembly454largecontigs.fna  asmion.largecontigs.stats q15. which of the technologies provided the best assembly of the e. coli strain? congrat lations yo finished the exercise!