Dinh/Dinh 2013/NOTES/2013-6-24

From ZhangLabWiki
Revision as of 19:49, 26 June 2013 by >Dinh
Jump to navigation Jump to search

Analysis of Blueprint MiSeq Test Run

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

Mandatory sites


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