Dinh/Dinh 2013/NOTES/2013-6-24: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
>Dinh
>Dinh
Line 280: Line 280:
* Calculate probe bias to re-subset probes as in Deng 2009 NBT paper.
* 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.
* Also, identify issues with the zero or low coverage mandatory and recommended sites and re-design probes.
===Re-design and try again with these sites===
===Identify problematic probes for mandatory and recommended targets===
* Note: recommended_29 and recommended_30 became optional.
* Note: recommended_29 and recommended_30 became optional.
* '''Re-design both strands for:'''
* '''Re-design both strands for:'''

Revision as of 01:49, 27 June 2013

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

  • 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.

Identify problematic probes for mandatory and recommended targets

  • Note: recommended_29 and recommended_30 became optional.
  • Re-design both strands for:
 *mandatory_1 (too in the middle?)
 mandatory_11 
 *mandatory_15 (too in the middle?)
 mandatory_16
 recommended_4 (did not designed probe for)
 recommended_21
 recommended_20 (low efficiency)
 recommended_27 (low efficiency)
  • Re-design Watson only
 recommended_2
 recommended_5 
 recommended_8
 recommended_22
 recommended_25
 mandatory_2 (low efficiency)
 mandatory_8 (low efficiency)
 mandatory_10 (low efficiency)
 mandatory_14 (low efficiency)
 recommended_9 (low efficiency)
 recommended_23 (low efficiency)
  • Re-design Crick only
 mandatory_5 
 recommended_16
 mandatory_4 (low efficiency)
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	+

order locus_identifier chrom probe_id probe_strand Crick total Watson total Capture
3 mandatory_3 chr4 cg14372037 + 1711 83 Both
6 mandatory_6 chr17 cg00960700 - 476 46 Both
7 mandatory_7 chr3 cg19442495 + 290 122 Both
9 mandatory_9 chr7 cg07945582 + 1346 486 Both
12 mandatory_12 chr2 cg12630082 + 107 242 Both
13 mandatory_13 chr15 cg07382129 + 194 4162 Both
17 recommended_1 chr5 cg06837426 - 584 43 Both
19 recommended_3 chr19 cg27321876 - 1345 102 Both
22 recommended_6 chr10 cg09060610 - 3789 61 Both
23 recommended_7 chr17 cg15298719 - 2017 77 Both
26 recommended_10 chr1 cg09846895 - 303 2362 Both
27 recommended_11 chr21 cg05931989 - 400 2966 Both
28 recommended_12 chr2 cg09590377 - 16 39 Both
29 recommended_13 chr20 cg23635789 + 448 1050 Both
30 recommended_14 chr17 cg03048083 - 4300 222 Both
31 recommended_15 chr16 cg26648465 - 387 1073 Both
33 recommended_17 chr17 cg05306745 + 69 337 Both
34 recommended_18 chr11 cg25169679 - 560 535 Both
35 recommended_19 chr16 cg07197480 - 4091 494 Both
40 recommended_24 chr1 cg20145149 - 4970 417 Both
42 recommended_26 chr1 cg07382920 - 320 500 Both
44 recommended_28 chr17 cg04819499 - 569 2377 Both
47 recommended_31 chrX cg24347720 + 5337 378 Both
48 recommended_32 chr17 cg21824343 + 2088 593 Both
49 recommended_33 chr8 cg00552087 - 263 15227 Both
50 recommended_34 chr1 cg24019564 - 11625 489 Both
5 mandatory_5 chr2 cg19426625 - 0 1695 Watson only
32 recommended_16 chr4 cg12311175 + 0 45 Watson only
36 recommended_20 chr8 cg17316966 + 8 0 Low both
43 recommended_27 chr16 cg01128460 - 1 9 Low both
2 mandatory_2 chr1 cg21646186 + 51 8 Low Watson
8 mandatory_8 chr7 cg06310157 - 1719 11 Low Watson
10 mandatory_10 chr7 cg15025536 + 3067 1 Low Watson
14 mandatory_14 chr4 cg01572513 - 509 11 Low Watson
25 recommended_9 chr13 cg22583065 + 913 6 Low Watson
39 recommended_23 chr6 cg14173662 - 1651 3 Low Watson
4 mandatory_4 chr2 cg04457196 + 3 30 Low Crick
1 mandatory_1 chr4 cg22178613 - 0 0 None
11 mandatory_11 chr7 cg12743416 - 0 0 None
15 mandatory_15 chr17 cg15044248 + 0 0 None
16 mandatory_16 chr7 cg24084681 + 0 0 None
20 recommended_4 chr16 cg06405517 - 0 0 None
37 recommended_21 chr10 cg00660608 - 0 0 None
45 recommended_29 chr4 cg22475974 + 0 0 None
46 recommended_30 chr7 cg18278817 - 0 0 None
18 recommended_2 chr20 cg00648301 + 134 0 Crick only
21 recommended_5 chr7 cg18680977 - 4863 0 Crick only
24 recommended_8 chr2 cg05986044 + 63 0 Crick only
38 recommended_22 chr7 cg09361748 - 698 0 Crick only
41 recommended_25 chr12 cg21898944 - 335 0 Crick only