Sam:LabNotes/Microbiome-new/2010-12-8: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
>Sam Chiang
>Sam Chiang
 
(17 intermediate revisions by the same user not shown)
Line 1: Line 1:
='''General statiscs (whole reads) and contamination analysis (failed to aligned reads) - HC-MDA-Ecoli-SAGs-Ind1to8'''=
='''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
Line 33: Line 33:


*Perform BLAST on filtered contigs
*Perform BLAST on filtered contigs
**BLAST setting: -K 1 -m 8 -a 2 -v 1 –b 1(using -m 8 output format to extract accession number)
**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)
**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   
***'''NOTE''': Use -K 1 -v 1 -b 1 for reporting only one hit   
**Taxonomic analysis by NCBI Entrez (against nucletide database)
**Taxonomic analysis by NCBI Batch Entrez (against nucletide database). Transform GI number to UID for taxonomic analysis


===Results===
===Results===
Line 86: Line 86:
*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.
*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.