Dinh/Dinh 2013/NOTES/2013-6-24: Difference between revisions
Jump to navigation
Jump to search
>Dinh m (→Probe Bias) |
>Dinh m (→Probe Bias) |
||
Line 278: | Line 278: | ||
|} | |} | ||
==Probe Bias== | ==Probe Bias== | ||
* Calculate probe bias to re-subset probes as in Deng 2009 NBT paper. | |||
* Also, identify issues with the zero or low coverage mandatory and recommended sites and re-design probes. | |||
===Mandatory sites=== | ===Mandatory sites=== | ||
<nowiki> | <nowiki> |
Revision as of 23:49, 26 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
kmerfreq -k 17 -t 4 -c -1 -p kmer_table read.lst >kmerfreq.cout 2>kmerfreq.cerr
- 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
- f1 = COPE combined PE reads, f2 = leftover R1, f3 = leftover R2.
- 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
- Modified Prof. Zhang's code to remove clonal reads using UMI.
Final Blueprint capture experiment results
- The capture setup for Index 6 provided the best results.
- Mapping was performed without adapter end-trimming.
- capture specificity = # on target bp / total usable bp
- capture sensitivity = # target bp at 1x / target size
- capture enrichment = # on target bp / target size
index | total PE reads | total bps | total bps mapped | % bps mapped | total bps after clonal removal | % bps clonal | genome bp | average genome depth of coverage | total CpG depth of coverage | number of CpGs (1x) | average CpGs depth of coverage | capture specificity | capture sensitivity | capture enrichment | f/r correlation | f/r correlation (w/o clonal removal) | #mandatory CpG (max 16) | #recommended CpG (max 32) | #optional CpG (max 1024) |
1 | 546,091 | 273,045,500 | 109,362,218 | 40.05% | 13,695,568 | 87.48% | 439,410 | 31 | 736,350 | 16,720 | 44 | 82.77% | 55.12% | 35.07 | 0.8204 | 0.7376 | 8 | 22 | 675 |
2 | 1,201,482 | 600,741,000 | 223,904,642 | 37.27% | 31,069,492 | 86.12% | 700,376 | 44 | 1,638,154 | 23,633 | 69 | 83.07% | 64.59% | 79.85 | 0.7707 | 0.7222 | 10 | 23 | 786 |
3 | 2,107,061 | 1,053,530,500 | 337,353,528 | 32.02% | 59,858,799 | 82.26% | 1,448,705 | 41 | 3,082,267 | 42,205 | 73 | 80.03% | 73.22% | 148.22 | 0.8222 | 0.7619 | 10 | 27 | 869 |
4 | 1,576,107 | 788,053,500 | 276,382,437 | 35.07% | 86,119,440 | 68.84% | 2,394,363 | 36 | 4,322,851 | 67,792 | 64 | 82.13% | 78.36% | 218.84 | 0.8762 | 0.8343 | 12 | 29 | 920 |
5 | 2,416,376 | 1,208,188,000 | 497,412,146 | 41.17% | 88,991,063 | 82.11% | 1,618,255 | 55 | 4,617,178 | 47,342 | 98 | 84.13% | 75.74% | 231.64 | 0.8639 | 0.8170 | 11 | 27 | 898 |
6 | 2,025,366 | 1,012,683,000 | 386,996,516 | 38.21% | 127,581,738 | 67.03% | 3,036,498 | 42 | 6,328,115 | 81,327 | 78 | 84.30% | 81.11% | 332.75 | 0.9089 | 0.8677 | 12 | 30 | 946 |
7 | 1,632,245 | 816,122,500 | 396,514,234 | 48.59% | 116,341,160 | 70.66% | 1,499,425 | 78 | 6,022,204 | 44,654 | 135 | 90.63% | 78.05% | 326.23 | 0.9029 | 0.8467 | 12 | 28 | 919 |
8 | 2,699,273 | 1,349,636,500 | 614,958,148 | 45.56% | 79,887,075 | 87.01% | 1,314,474 | 61 | 4,271,697 | 38,539 | 111 | 87.39% | 72.11% | 216.01 | 0.8372 | 0.7962 | 11 | 27 | 852 |
Probe Bias
- Calculate probe bias to re-subset probes as in Deng 2009 NBT paper.
- Also, identify issues with the zero or low coverage mandatory and recommended sites and re-design probes.
Mandatory sites
reads site str 23892 mandatory_13:chr15:100249085-100249289 + 12876 mandatory_10:chr7:3025508-3025737 - 5540 mandatory_3:chr4:7526550-7526760 - 4183 mandatory_15:chr17:75369137-75369338 - 3944 mandatory_1:chr4:154710370-154710640 + 2965 mandatory_8:chr7:140218006-140218261 - 2579 mandatory_5:chr2:9518172-9518402 + 1839 mandatory_6:chr17:80709116-80709327 - 1358 mandatory_9:chr7:26206450-26206690 - 805 mandatory_14:chr4:147557728-147557999 - 633 mandatory_9:chr7:26206446-26206652 + 561 mandatory_12:chr2:42275654-42275925 + 537 mandatory_7:chr3:142837896-142838117 - 484 mandatory_7:chr3:142837897-142838117 + 331 mandatory_12:chr2:42275614-42275844 - 231 mandatory_13:chr15:100249082-100249302 - 230 mandatory_3:chr4:7526557-7526773 + 183 mandatory_6:chr17:80709042-80709272 + 67 mandatory_2:chr1:110052333-110052595 - 32 mandatory_4:chr2:164593111-164593326 + 17 mandatory_4:chr2:164593173-164593378 - 16 mandatory_14:chr4:147557730-147558008 + 9 mandatory_8:chr7:140218000-140218270 + 5 mandatory_2:chr1:110052282-110052542 + 5 mandatory_1:chr4:154710373-154710633 - 5 mandatory_10:chr7:3025518-3025736 + 0 mandatory_15:chr17:75369134-75369334 + 0 mandatory_5:chr2:9518220-9518443 - 0 mandatory_11:chr7:138229771-138230020 + 0 mandatory_11:chr7:138229843-138230113 - 0 mandatory_16:chr7:93520143-93520413 -
Recommended
reads site str 65845 recommended_33:chr8:3316771-3316991 + 19194 recommended_34:chr1:25257453-25257723 - 12188 recommended_14:chr17:43339401-43339611 - 12048 recommended_11:chr21:47783987-47784263 + 11668 recommended_24:chr1:115124392-115124598 - 11231 recommended_31:chrX:135333528-135333788 - 11132 recommended_5:chr7:27154944-27155156 - 10747 recommended_19:chr16:1017670-1017944 - 9879 recommended_6:chr10:7451206-7451436 - 7055 recommended_32:chr17:1633600-1633820 - 6674 recommended_10:chr1:161442627-161442877 + 5537 recommended_28:chr17:43044892-43045132 + 4945 recommended_15:chr16:85676275-85676545 + 4350 recommended_23:chr6:32042899-32043159 - 4311 recommended_3:chr19:38085548-38085818 - 2204 recommended_9:chr13:103052745-103053015 - 2117 recommended_13:chr20:36013306-36013546 + 2010 recommended_32:chr17:1633600-1633820 + 1950 recommended_22:chr7:71682175-71682435 - 1822 recommended_1:chr5:54281139-54281408 - 1516 recommended_19:chr16:1017669-1017943 + 1285 recommended_7:chr17:7165161-7165397 - 1185 recommended_13:chr20:36013306-36013536 - 1163 recommended_18:chr11:69197082-69197351 - 981 recommended_11:chr21:47783999-47784254 - 934 recommended_33:chr8:3316763-3316986 - 900 recommended_26:chr1:3567532-3567752 - 773 recommended_28:chr17:43044890-43045150 - 720 recommended_17:chr17:76921742-76921982 + 714 recommended_26:chr1:3567513-3567773 + 676 recommended_10:chr1:161442625-161442877 - 627 recommended_25:chr12:130589110-130589320 - 504 recommended_24:chr1:115124345-115124605 + 448 recommended_18:chr11:69197079-69197355 + 410 recommended_7:chr17:7165139-7165351 + 354 recommended_3:chr19:38085549-38085819 + 297 recommended_2:chr20:20349143-20349401 - 241 recommended_24:chr1:115124345-115124605 - 206 recommended_17:chr17:76921738-76921961 - 123 recommended_15:chr16:85676279-85676544 - 119 recommended_12:chr2:8597284-8597533 + 94 recommended_6:chr10:7451190-7451465 + 91 recommended_8:chr2:9614480-9614708 - 88 recommended_1:chr5:54281139-54281409 + 58 recommended_34:chr1:25257451-25257727 + 53 recommended_31:chrX:135333568-135333778 + 43 recommended_16:chr4:118975389-118975609 + 39 recommended_12:chr2:8597159-8597438 - 16 recommended_27:chr16:1243456-1243666 + 8 recommended_9:chr13:103052712-103052935 + 7 recommended_20:chr8:74878400-74878650 - 4 recommended_27:chr16:1243426-1243663 - 2 recommended_23:chr6:32042893-32043163 + 0 recommended_21:chr10:125853024-125853277 - 0 recommended_25:chr12:130588995-130589226 + 0 recommended_32:chr17:1633599-1633825 + 0 recommended_14:chr17:43339405-43339665 + 0 recommended_2:chr20:20349125-20349399 + 0 recommended_8:chr2:9614470-9614722 + 0 recommended_16:chr4:118975389-118975619 - 0 recommended_5:chr7:27154941-27155143 + 0 recommended_22:chr7:71682176-71682426 + 0 recommended_20:chr8:74878348-74878592 +