Sam:LabNotes/Microbiome-new/2010-12-8: Difference between revisions
Jump to navigation
Jump to search
>Sam Chiang |
>Sam Chiang |
||
(42 intermediate revisions by the same user not shown) | |||
Line 1: | Line 1: | ||
=''' | ='''General statiscs (whole reads) and contamination analysis (failed-to-align reads) - HC-MDA-Ecoli-SAGs-Ind1to8'''= | ||
==Background== | ==Background== | ||
Line 12: | Line 12: | ||
*Whole reads alignment using Bowtie => genome coverage calculate | *Whole reads alignment using Bowtie => genome coverage calculate | ||
**Script: Read2coverage-kz-3 | **Script: Read2coverage-kz-3 | ||
**Data resource: HC-MDA-Ecoli-SAGs-Ind1to8 | **Data resource: HC-MDA-Ecoli-SAGs-Ind1to8 (SE reads data set) | ||
**Template:E.coli-k12-mg1655 template | **Template:E.coli-k12-mg1655 template | ||
**Bowtie setting: -k 1 -l 28 -n 2 -m 1 --best --best --phred64-quals | **Bowtie setting: -k 1 -l 28 -n 2 -m 1 --best --best --phred64-quals | ||
**Downstream: Data was used for genome coverage calculation. | **Downstream: Data was used for genome coverage calculation. | ||
[[File:sam120610-wholereadbotie.jpg|800px]] | [[File:sam120610-wholereadbotie-2.jpg|800px]] | ||
*Using the same alignment processing, I extract the failed-to-aligned and non-specific reads into different folder. | *Using the same alignment processing, I extract the failed-to-aligned and non-specific reads into different folder. | ||
Line 26: | Line 26: | ||
*ABySS de novo assembly of failed-to-align reads with (using single-end assembly with k-mer 27) | *ABySS de novo assembly of failed-to-align reads with (using single-end assembly with k-mer 27) | ||
**Using abyss-fac function to calculate statistics | **Using abyss-fac function to calculate statistics | ||
'''ABySS-fac report''' | |||
[[File:sam120610-abyss assembly of fail-to-aligne reads.jpg|700px]] | [[File:sam120610-abyss assembly of fail-to-aligne reads.jpg|700px]] | ||
* | *Length filtering and keep only >300bp contigs. | ||
**Script: seqLenFilger.pl | **Script: seqLenFilger.pl | ||
*Perform BLAST on filtered contigs | *Perform BLAST on filtered contigs | ||
**BLAST setting: -K 1 -m 8 -a 2 (using -m 8 output format to extract | **BLAST setting: -K 1 -m 8 -a 2 -v 1 –b 1(using -m 8 output format to extract GI number) | ||
**BLAST setting: -K 1 -m 8 -a 2 (using -m 9 output format to get BLAST statistics) | **BLAST setting: -K 1 -m 8 -a 2 -v 1 –b 1(using -m 9 output format to get BLAST statistics) | ||
**Taxonomic analysis by NCBI | ***'''NOTE''': Use -K 1 -v 1 -b 1 for reporting only one hit | ||
**Taxonomic analysis by NCBI Batch Entrez (against nucletide database). Transform GI number to UID for taxonomic analysis | |||
===Results=== | ===Results=== | ||
Line 47: | Line 50: | ||
'''E.coli taxonomic ranks''' | |||
------------------------- | |||
Domain: Eubacteria; | |||
'''E.coli | |||
--------------------- | |||
Domain: | |||
Phylum: Proteobacteria; | Phylum: Proteobacteria; | ||
Class: Gammaproteobacteria; | Class: Gammaproteobacteria; | ||
Order: Enterobacteriales; | Order: Enterobacteriales; | ||
Family: Enterobacteriaceae; | Family: Enterobacteriaceae; | ||
Genus: Escherichia | Genus: Escherichia; | ||
Species: E. coli | Species: E. coli | ||
=== | '''Yeast taxonomic ranks''' - (to detect if Yeast DNA contamination was introduced through Yeast tRNA carrier) | ||
------------------------ | |||
Super kindom: Eukaryota; | |||
Domain: Fungi; | |||
Phylum: Ascomycota; | |||
Class: Saccharomycetes; | |||
Order: Saccharomycetales; | |||
Family: Saccharomycetaceae; | |||
Genus: Saccharomyces; | |||
Species: Saccharomyces cerevisiae | |||
===Discussions=== | |||
*Most of fail-to-align reads can only be assembled to short contigs (under k-mer 27) | *Most of fail-to-align reads can only be assembled to short contigs (under k-mer 27) | ||
*The contigs assembled from fail-to-aligned reads(against E.coli template) contain a significant amount of human DNA contigs. | *The contigs assembled from fail-to-aligned reads(against E.coli template) '''contain a significant amount of human DNA contigs'''. | ||
* | ID1 ID2 ID3 ID4 ID5 ID6 ID7 ID8 | ||
**Sequencing error => mis-match >3bp => fail to meet mapping criteria in bowtie? | ----------------------------------------------------------- | ||
**Chimeric junction in these "single-end" reads? | Homo hits N/A 22 0 25 155 86 35 83 | ||
All hits N/A 75 15 106 200 145 56 120 | |||
(%) 29.3% 23.6% 77.5% 59.3% 62.5% 69.1% | |||
All failed-to-align contig hit | |||
*'''A certain level fail-to-align contigs also match E.coli genome (e.g. ID-3 and ID-8)'''. It is unknown why some E.coli reads are failed to map to E.coli genome in bowtie. | |||
***Sequencing error => mis-match >3bp => fail to meet mapping criteria in bowtie? | |||
***Chimeric junction in these "single-end" reads? | |||
*ID-1 showed no BLAST hit using length>300 contigs. I tried using length>200 contigs for BLAST but it still came out no hit. | *ID-1 showed no BLAST hit using length>300 contigs. I tried using length>200 contigs for BLAST but it still came out no hit. | ||
*I am thinking that maybe it's a better idea to '''use MEGAN metagenomics tool for analyzing these fail-to-align reads''' directly, which it can draw a taxonomic distribution plot easily. | |||
==Testing using less-stringent report suppression setting on bowtie mapping== | |||
*I noticed that the genome coverage calculation script (Read2coverage-kz-3) calculate the coverage from bowtie.out file, however the bowtie.out file only report unique mapped reads (confirmed on http://genome-tech.ucsd.edu/LabNotes/index.php/Sam:LabNotes/Microbiome-new/2010-12-13). | |||
**I guess there are genome coverages not counted if they are from non-unique mapped reads, since non-uniqe mapped reads were not reported in bowtie.out. | |||
===Re-analyze genome coverage using new scripts with less-stringent report suppression=== | |||
*test data set: s_4_ID6.txt from SE data set (HL083 E.coli SAGs) | |||
**This data set has the highest genome coverage than others. | |||
*'''Read2coverage-kz-3-suppress_5''' | |||
**(allow reporting 5 alignments from 1 read) | |||
*Results | |||
# reads processed: 2372473 | |||
# reads with at least one reported alignment: 1512449 (63.75%) | |||
# reads that failed to align: 512805 (21.61%) | |||
# reads with alignments suppressed due to -m: 347219 (14.64%) | |||
Reported 1512449 alignments to 1 output stream(s) | |||
total bps coverred = 1824363 | |||
*'''Read2coverage-kz-3-suppress_10''' | |||
**(allow reporting 10 alignments from 1 read) | |||
*Results | |||
# reads processed: 2372473 | |||
# reads with at least one reported alignment: 1857421 (78.29%) | |||
# reads that failed to align: 512805 (21.61%) | |||
# reads with alignments suppressed due to -m: 2247 (0.09%) | |||
Reported 1857421 alignments to 1 output stream(s) | |||
total bps coverred = 1869404 | |||
==Conclusions== | |||
*Test data: E.coli SAG HL083(SE) s_4_ID6.txt | |||
{| {{table}} | |||
| align="center" style="background:#f0f0f0;"|'''Report suppressor''' | |||
| align="center" style="background:#f0f0f0;"|'''m = 1''' | |||
| align="center" style="background:#f0f0f0;"|'''m = 5''' | |||
| align="center" style="background:#f0f0f0;"|''' m= 10''' | |||
|- | |||
| proccessed reads||2372473||2372473||2372473 | |||
|- | |||
| unique alignment||1411597 (59.50%)||1512449 (63.75%)||1857421 (78.29%) | |||
|- | |||
| failed-to-align||512805 (21.61%)||512805 (21.61%)||512805 (21.61%) | |||
|- | |||
| non-unique alignment||448071 (18.89%)||347219 (14.64%)||2247 (0.09%) | |||
|- | |||
| Total covered bp||1778296||1824363||1869404 | |||
|- | |||
| Coverage (E.coli K12)||38.33%||39.32%||40.30% | |||
|- | |||
| | |||
|} | |||
*Tune down the report suppression strigency did work on increasing genome coverage, but the effect is limited. Even with m=10 setting where alomst all non-unique alignments (18%) were reported (counted in genome coverage), the genome coverage was only increased by ~2% at most. | |||
*'''The genome coverage is mainly contributed by unique alignment'''. | |||
==Identify where did those non-unique reads map== | |||
*Perform ABYSS on non-unique read of s_4_ID6.nonspecific reads | |||
*I mapped the assembled contigs on E.coli reference genome using BLAT, and found that these contigs are randomly distributed. |
Latest revision as of 01:22, 12 August 2011
General statiscs (whole reads) and contamination analysis (failed-to-align reads) - HC-MDA-Ecoli-SAGs-Ind1to8[edit]
Background[edit]
- Samples were prepared from -(short and long amplified E.coli SAGs -3rd test 11-07-'10)
- Eight E.coli SAGs samples were made to Illumina sequencing libraries using (N2)Primer-based barcodes ID1 to ID8.
- These libraries were equally pulled (20ng per library) and sequenced in pair-end(PE) mode.
- The sequenced data processed by Dr. Zhang for base-calling and saved into FASTQ files.
- Since I don't have genome-tech server permission, I asked Dinh's help to transfer data into Genemapster server.
Analysis approaches[edit]
Single-end data analysis[edit]
- Whole reads alignment using Bowtie => genome coverage calculate
- Script: Read2coverage-kz-3
- Data resource: HC-MDA-Ecoli-SAGs-Ind1to8 (SE reads data set)
- Template:E.coli-k12-mg1655 template
- Bowtie setting: -k 1 -l 28 -n 2 -m 1 --best --best --phred64-quals
- Downstream: Data was used for genome coverage calculation.
File:Sam120610-wholereadbotie-2.jpg
- Using the same alignment processing, I extract the failed-to-aligned and non-specific reads into different folder.
- Script: Read2coverage-kz4
- Data resource: HC-MDA-Ecoli-SAGs-Ind1to8
- Template:E.coli-k12-mg1655 template
- Bowtie setting: -k 1 -l 28 -n 2 -m 1 --best --best --phred64-quals --un [failed-to-aligned file output] --max [nonspecific file output]
- ABySS de novo assembly of failed-to-align reads with (using single-end assembly with k-mer 27)
- Using abyss-fac function to calculate statistics
ABySS-fac report File:Sam120610-abyss assembly of fail-to-aligne reads.jpg
- Length filtering and keep only >300bp contigs.
- Script: seqLenFilger.pl
- Perform BLAST on filtered contigs
- BLAST setting: -K 1 -m 8 -a 2 -v 1 –b 1(using -m 8 output format to extract GI number)
- BLAST setting: -K 1 -m 8 -a 2 -v 1 –b 1(using -m 9 output format to get BLAST statistics)
- NOTE: Use -K 1 -v 1 -b 1 for reporting only one hit
- Taxonomic analysis by NCBI Batch Entrez (against nucletide database). Transform GI number to UID for taxonomic analysis
Results[edit]
px300 px300 Details of blast hits (PDF) Media:ID2-failed-to-align-L300-accession.pdf Media:ID3-failed-to-align-L300-accession.pdf Media:ID4-failed-to-align-L300-accession.pdf Media:ID5-failed-to-align-L300-accession.pdf Media:ID6-failed-to-align-L300-accession.pdf Media:ID7-failed-to-align-L300-accession.pdf Media:ID8-failed-to-align-L300-accession.pdf
E.coli taxonomic ranks ------------------------- Domain: Eubacteria; Phylum: Proteobacteria; Class: Gammaproteobacteria; Order: Enterobacteriales; Family: Enterobacteriaceae; Genus: Escherichia; Species: E. coli
Yeast taxonomic ranks - (to detect if Yeast DNA contamination was introduced through Yeast tRNA carrier) ------------------------ Super kindom: Eukaryota; Domain: Fungi; Phylum: Ascomycota; Class: Saccharomycetes; Order: Saccharomycetales; Family: Saccharomycetaceae; Genus: Saccharomyces; Species: Saccharomyces cerevisiae
Discussions[edit]
- Most of fail-to-align reads can only be assembled to short contigs (under k-mer 27)
- The contigs assembled from fail-to-aligned reads(against E.coli template) contain a significant amount of human DNA contigs.
ID1 ID2 ID3 ID4 ID5 ID6 ID7 ID8 ----------------------------------------------------------- Homo hits N/A 22 0 25 155 86 35 83 All hits N/A 75 15 106 200 145 56 120 (%) 29.3% 23.6% 77.5% 59.3% 62.5% 69.1%
All failed-to-align contig hit
- A certain level fail-to-align contigs also match E.coli genome (e.g. ID-3 and ID-8). It is unknown why some E.coli reads are failed to map to E.coli genome in bowtie.
- Sequencing error => mis-match >3bp => fail to meet mapping criteria in bowtie?
- Chimeric junction in these "single-end" reads?
- ID-1 showed no BLAST hit using length>300 contigs. I tried using length>200 contigs for BLAST but it still came out no hit.
- I am thinking that maybe it's a better idea to use MEGAN metagenomics tool for analyzing these fail-to-align reads directly, which it can draw a taxonomic distribution plot easily.
Testing using less-stringent report suppression setting on bowtie mapping[edit]
- I noticed that the genome coverage calculation script (Read2coverage-kz-3) calculate the coverage from bowtie.out file, however the bowtie.out file only report unique mapped reads (confirmed on http://genome-tech.ucsd.edu/LabNotes/index.php/Sam:LabNotes/Microbiome-new/2010-12-13).
- I guess there are genome coverages not counted if they are from non-unique mapped reads, since non-uniqe mapped reads were not reported in bowtie.out.
Re-analyze genome coverage using new scripts with less-stringent report suppression[edit]
- test data set: s_4_ID6.txt from SE data set (HL083 E.coli SAGs)
- This data set has the highest genome coverage than others.
- Read2coverage-kz-3-suppress_5
- (allow reporting 5 alignments from 1 read)
- Results
# reads processed: 2372473 # reads with at least one reported alignment: 1512449 (63.75%) # reads that failed to align: 512805 (21.61%) # reads with alignments suppressed due to -m: 347219 (14.64%) Reported 1512449 alignments to 1 output stream(s) total bps coverred = 1824363
- Read2coverage-kz-3-suppress_10
- (allow reporting 10 alignments from 1 read)
- Results
# reads processed: 2372473 # reads with at least one reported alignment: 1857421 (78.29%) # reads that failed to align: 512805 (21.61%) # reads with alignments suppressed due to -m: 2247 (0.09%) Reported 1857421 alignments to 1 output stream(s) total bps coverred = 1869404
Conclusions[edit]
- Test data: E.coli SAG HL083(SE) s_4_ID6.txt
Report suppressor | m = 1 | m = 5 | m= 10 |
proccessed reads | 2372473 | 2372473 | 2372473 |
unique alignment | 1411597 (59.50%) | 1512449 (63.75%) | 1857421 (78.29%) |
failed-to-align | 512805 (21.61%) | 512805 (21.61%) | 512805 (21.61%) |
non-unique alignment | 448071 (18.89%) | 347219 (14.64%) | 2247 (0.09%) |
Total covered bp | 1778296 | 1824363 | 1869404 |
Coverage (E.coli K12) | 38.33% | 39.32% | 40.30% |
- Tune down the report suppression strigency did work on increasing genome coverage, but the effect is limited. Even with m=10 setting where alomst all non-unique alignments (18%) were reported (counted in genome coverage), the genome coverage was only increased by ~2% at most.
- The genome coverage is mainly contributed by unique alignment.
Identify where did those non-unique reads map[edit]
- Perform ABYSS on non-unique read of s_4_ID6.nonspecific reads
- I mapped the assembled contigs on E.coli reference genome using BLAT, and found that these contigs are randomly distributed.