Dinh/Dinh 2013/NOTES/2013-6-24: Difference between revisions
Jump to navigation
Jump to search
>Dinh mNo edit summary |
>Dinh mNo edit summary |
||
Line 168: | Line 168: | ||
===Mapping statistics=== | ===Mapping statistics=== | ||
* With clipping, we got 2.846 Gbps mapped, without clipping, we got 2.842 Gbps mapped. Clipping did not improve the amount of usable bp by a lot. | * With clipping, we got 2.846 Gbps mapped, without clipping, we got 2.842 Gbps mapped. Clipping did not improve the amount of usable bp by a lot. | ||
* Libraries Index 7 and Index 8 have the best mapping rates. | |||
{| {{table}} | {| {{table}} | ||
| align="center" style="background:#f0f0f0;"|'''index''' | | align="center" style="background:#f0f0f0;"|'''index''' |
Revision as of 19:26, 25 June 2013
Analysis of Blueprint MiSeq Test Run
- 2 X 250bp MiSeq PE run, total = 13,117,847 PE reads
- Noi's labnote page: http://genome-tech.ucsd.edu/LabNotes/index.php/Noi/NOTES/2013-6-6
Pre-processing
Standard trimming + UMI
- Obtain UMI from the first 10bp of read 1, label both reads with UMI
- Trim 27 bp from 5 prime end of read 1 and read 2
Adapter removal with fastq-mcf
- Remove adapter sequences using fastq-mcf ?
- Since the insert = (target_size= 200-280) + (2 arms ~ 64 bp) + (UMI = 10 bp) = 274-354bp, we would not have sequenced the adapters with just 250bp from each end.
- Try removing adapters because mapping rate was ~40%.
- Adapters:
>Linker TTGGAGGCTCATCGTTCCTATTCAGGCAGATGTTATCGAGGTCCGAC >Linker_rev GTCGGACCTCGATAACATCTGCCTGAATAGGAACGATGAGCCTCCAA
fastq-mcf results
- Note that only R1 was end-trimmed using "Linker" and only R2 was end-trimmed using "Linker_rev"
index | Total reads | Clipped 'end' reads | % clipped | Too short after clip | % too short |
1, R1 | 546,091 | 265,765 | 48.67% | 16,468 | 3.02% |
1, R2 | 546,091 | 330,822 | 60.58% | 16,112 | 2.95% |
2, R1 | 1,201,482 | 599,205 | 49.87% | 46,756 | 3.89% |
2, R2 | 1,201,482 | 733,958 | 61.09% | 45,880 | 3.82% |
3, R1 | 2,107,061 | 1,159,863 | 55.05% | 30,221 | 1.43% |
3, R2 | 2,107,061 | 1,348,459 | 64.00% | 29,076 | 1.38% |
4, R1 | 1,576,107 | 821,955 | 52.15% | 16,478 | 1.05% |
4, R2 | 1,576,107 | 1,000,442 | 63.48% | 16,703 | 1.06% |
5, R1 | 2,416,376 | 1,172,370 | 48.52% | 30,785 | 1.27% |
5, R2 | 2,416,376 | 1,470,535 | 60.86% | 29,880 | 1.24% |
6, R1 | 2,025,366 | 1,044,010 | 51.55% | 17,056 | 0.84% |
6, R2 | 2,025,366 | 1,266,104 | 62.51% | 16,306 | 0.81% |
7, R1 | 1,632,245 | 737,076 | 45.16% | 12,924 | 0.79% |
7, R2 | 1,632,245 | 398,506 | 24.41% | 12,788 | 0.78% |
8, R1 | 2,699,273 | 1,252,029 | 46.38% | 24,410 | 0.90% |
8, R2 | 2,699,273 | 668,800 | 24.78% | 22,635 | 0.84% |
Merge R1 and R2 with COPE
- Create kmer_table for COPE
- Use COPE to combine overlapping read 1 and read 2
for f in 1Index1_S1 1Index2_S2 1Index3_S3 2Index4_S6 1Index5_S4 2Index6_S7 2Index7_S8 1Index8_S5 do ./getUMI.pl 130620_MiSeq/TES1-${f}_L001_R1_001.fastq 130620_MiSeq/TES1-${f}_L001_R2_001.fastq $f ~/softwares/cope-src-v1.1.3/src/cope -a $f.R1.fq -b $f.R2.fq -o $f.fq -2 $f.leftR1.fq -3 $f.leftR2.fq -m 1 -t kmer_table.freq.cz -f kmer_table.freq.cz.len >cope.$f.log 2>cope.$f.error rm $f.R1.fq $f.R2.fq done
COPE results
- Without adapter removal:
index | total_pairs | connected_pairs | connect_ratio(%) | low_quality_pairs | low_quality_ratio(%) |
1 | 546,091 | 37,624 | 6.88969 | 239413 | 43.8412 |
2 | 1,201,482 | 76,332 | 6.35315 | 603335 | 50.2159 |
3 | 2,107,061 | 120,889 | 5.73733 | 1155859 | 54.8565 |
5 | 2,416,376 | 157,448 | 6.51587 | 1124046 | 46.5178 |
8 | 2,699,273 | 191,332 | 7.08828 | 1138412 | 42.1748 |
4 | 1,576,107 | 100,314 | 6.36467 | 805229 | 51.0897 |
6 | 2,025,366 | 131,396 | 6.48752 | 990421 | 48.9008 |
7 | 546,091 | 35,800 | 6.55568 | 239459 | 43.8497 |
- With adapter (linker sequence) removal:
index | total_pairs | connected_pairs | connect_ratio(%) | low_quality_pairs | low_quality_ratio(%) |
1 | 529623 | 20770 | 3.92166 | 195892 | 36.9871 |
2 | 1154726 | 42045 | 3.64112 | 477352 | 41.339 |
3 | 2076840 | 65820 | 3.16924 | 970053 | 46.7081 |
5 | 2385591 | 96387 | 4.04038 | 964834 | 40.4442 |
8 | 2674863 | 126170 | 4.71688 | 1012342 | 37.8465 |
4 | 1559629 | 55184 | 3.53828 | 684665 | 43.8992 |
6 | 2008310 | 76903 | 3.82924 | 820560 | 40.8582 |
7 | 1619321 | 80497 | 4.97103 | 584325 | 36.0846 |
Alignment with Bowtie2
- The long reads were not compatible with our pipeline using Bowtie2, because the aligner suppresses read name lines which are longer than 256 characters. This causes truncated original reads stored in the read name line.
- Added code to split long reads into multiple short reads to use Bowtie2
- File:SmartBisReadMapper.txt
- Mapping pipeline:
cur_dir="/media/3TB_Dinh/Blueprint" reads_dir="/media/3TB_Dinh/Blueprint" email="diep.hue.dinh@gmail.com" bisReadMapper="/home/ddiep/scripts/MethylationPipeline/scripts/smartBisReadMapper.pl" template_fwd="/media/2TB_storeA/BisRef/bisHg19/hg19.fa.bis.fwd.bowtie2" template_rev="/media/2TB_storeA/BisRef/bisHg19/hg19.fa.bis.rev.bowtie2" template_fa="/media/2TB_storeA/BisRef/bisHg19/hg19.fa" soap="/home/ddiep/softwares/soap2.21release/soap" bowtie="bowtie2" cpg="/media/2TB_storeA/BisRef/bisHg19/C_Pos/hg19.fa.cpg.positions.txt" INDX="1Index1_S1 1Index2_S2 1Index3_S3 2Index4_S6 1Index5_S4 2Index6_S7 2Index7_S8 1Index8_S5" # make sure the qual_base variable is set correctly cd $cur_dir for s in ${INDX} do f="$s.fq" g="$s.leftR1.fq" h="$s.leftR2.fq" n="$s-MS" mkdir $n echo "cd $cur_dir/$n" > $n.job echo "$bisReadMapper -r $reads_dir/$f -W $template_fwd -C $template_rev -g $template_fa -a $bowtie -p 16 -b 33 -n $n.f1 -q 20 -l $cpg > $n.f1.statusMbias 2>$n.f1.err" >> $n.job echo "$bisReadMapper -r $reads_dir/$g -W $template_fwd -C $template_rev -g $template_fa -a $bowtie -p 16 -b 33 -n $n.f2 -q 20 -l $cpg > $n.f2.statusMbias 2>$n.f2.err" >> $n.job echo "$bisReadMapper -r $reads_dir/$h -W $template_fwd -C $template_rev -g $template_fa -a $bowtie -p 16 -b 33 -n $n.f3 -q 20 -l $cpg > $n.f3.statusMbias 2>$n.f3.err" >> $n.job echo "rm *encoded" >> $n.job done
Mapping statistics
- With clipping, we got 2.846 Gbps mapped, without clipping, we got 2.842 Gbps mapped. Clipping did not improve the amount of usable bp by a lot.
- Libraries Index 7 and Index 8 have the best mapping rates.
index | read file | total clipped bases | total clipped mapped bases | %mapped(clipped) | total bases | total mapped bases | %mapped |
1 | f1 | 8255745 | 5487013 | 66.46% | 16078514 | 9619049 | 59.83% |
1 | f2 | 79964964 | 53864094 | 67.36% | 104318024 | 50945748 | 48.84% |
1 | f3 | 89749717 | 50800000 | 56.60% | 99998101 | 48797421 | 48.80% |
2 | f1 | 16566696 | 11031637 | 66.59% | 32606092 | 19771959 | 60.64% |
2 | f2 | 171105395 | 112357034 | 65.67% | 230586489 | 106229965 | 46.07% |
2 | f3 | 189062773 | 102110571 | 54.01% | 215612917 | 97902718 | 45.41% |
3 | f1 | 25104735 | 15365281 | 61.20% | 51734457 | 30631757 | 59.21% |
3 | f2 | 285378749 | 171910645 | 60.24% | 405686018 | 161655992 | 39.85% |
3 | f3 | 330732139 | 151852920 | 45.91% | 372674273 | 145065779 | 38.93% |
4 | f1 | 21403635 | 13378820 | 62.51% | 42833085 | 24714027 | 57.70% |
4 | f2 | 222330667 | 138200084 | 62.16% | 300568658 | 131048023 | 43.60% |
4 | f3 | 250430069 | 125847956 | 50.25% | 281078201 | 120620387 | 42.91% |
5 | f1 | 38452894 | 26019949 | 67.67% | 67355222 | 42886852 | 63.67% |
5 | f2 | 360645974 | 247314571 | 68.58% | 463290839 | 235545215 | 50.84% |
5 | f3 | 401191540 | 226901252 | 56.56% | 438788886 | 218980079 | 49.91% |
6 | f1 | 30067978 | 18916755 | 62.91% | 56145035 | 34232803 | 60.97% |
6 | f2 | 289462431 | 191584810 | 66.19% | 389224257 | 183501240 | 47.15% |
6 | f3 | 328112598 | 176045539 | 53.65% | 364481484 | 169262473 | 46.44% |
7 | f1 | 32538699 | 22178956 | 68.16% | 49348430 | 34280682 | 69.47% |
7 | f2 | 253380410 | 192127161 | 75.83% | 313644714 | 186792389 | 59.56% |
7 | f3 | 279659004 | 180487623 | 64.54% | 300856044 | 175441163 | 58.31% |
8 | f1 | 50524286 | 33292157 | 65.89% | 81821672 | 53601105 | 65.51% |
8 | f2 | 415435035 | 301628773 | 72.61% | 516798645 | 291900532 | 56.48% |
8 | f3 | 463261828 | 277764207 | 59.96% | 496294120 | 269456511 | 54.29% |
Remove clonal reads
- Use Prof. Shang's code to remove clonal reads using UMI