Dinh 2011/NOTES/2011-10-29
Jump to navigation
Jump to search
Specificity of DMR330K capture[edit]
- The latest best PE110 data is from Dr. Kang Zhang's group.
- Mapped one data set (86.39%) mapping rate.
- Ran bisPileup2coverage.pl :
########################### Total mappable bases = 2340766083 Total on-target bases = 2332481117 Specificity = 0.9965 ########################### total capture targets = 35829723 bp all_mapped on_target %target_covered DP>=1x 32120729 27091436 0.756116255768988 DP>=10x 20497181 19595293 0.546900488178488 DP>=20x 16075410 15554106 0.434111812698078
Further testing
- I've recently mapped a sample of 2.5 million reads from sample KEBR007A (UPenn48 dataset). Pileup files are on genome-miner: /home/ddiep/LibraryFree_SuppFig_Data/KEBR007A_1.txt.*.pileup
- I ran bisPileupPair2coverage.pl and got:
Total mappable bases = 77956846 Total on-target bases = 76934722 Specificity = 0.9869
- I modified bisPileupPair2coverage.pl to print out all positions counted as on-target (marked "rev:" for reverse strand positions): File:BisPileupPair2coverage DD2.txt
./bisPileupPair2coverage_DD2.pl CpG330KProbesInfo_nonOverlappingBED.txt KEBR007A_1.txt.fwd.pileup KEBR007A_1.txt.rev.pileup > tmp_ontarget awk '{if($0 ~ /chr/) print $0;}' tmp_ontarget | wc -l returned 8787716 awk '{if($0 ~ /chr/) print $0;}' tmp_ontarget | sort -u | wc -l returned 8787581 awk '{if($0 ~ /chr/) print $0;}' tmp_ontarget | sort | uniq -D > double_counted_targets
- There were some small number of duplicate counting which was not even by counting twice on Watson and Crick
ie. chr16 65518303 A A 30 0 30 1 ^?. g chr16 65518303 A A 30 0 30 1 ^?. g
- Found a mistake in CpG330KProbesInfo_nonOverlappingBED.txt file:
chr16 65518122 65518303 chr16 65518303 65518484 *The position 65518303 is covered twice.
- Fixed
Old non overlapping CpG330K bed file: 142747 regions Fixed non overlapping CpG330K bed file: 140560 regions
- Re-ran analysis:
./bisPileupPair2coverage_DD2.pl CpG330KProbesInfo_nonOverlappingBED-10292011.txt KEBR007A_1.txt.fwd.pileup KEBR007A_1.txt.rev.pileup > tmp_ontarget awk '{if($0 ~ /chr/) print $0;}' tmp_ontarget | sort | uniq -D > double_counted_targets wc -l double_counted_targets returns 0
Total mappable bases = 77956846 Total on-target bases = 76934722 Specificity = 0.9869
- Wrote simple script to count all on-target bases:
my $sum = 0; while(my $line = <STDIN>){ my @f =split /\t/, $line; next if(!$f[7]); $sum+=$f[7]; } print $sum,"\n";
perl countPileupBases.pl < tmp_ontarget returns 72431881
- Need to get the correct on-target base counts.
- Simple counting script: File:BisPileupPair2specificity.txt
./bisPileupPair2specificity.pl CpG330KProbesInfo_nonOverlappingBED-10292011.txt KEBR007A_1.txt.fwd.pileup KEBR007A_1.txt.rev.pileup OUTPUT: Size of target = 35827535 Total mappable bases = 77956846 Total target bases = 72431881 (0.929)
Fixed non overlapping LC4K bed file : 803 regions * Dr. Zhang's NM_LC_Probes_Filtered_09232011_non-overlapped.BED.txt have 804 regions. Comparing the two, there three entries were unique between the two files: In mine: chr16,55259870,55260287 In KZ: chr16,55259888,55260214 In KZ: chr16,55261003,55261181
The entries in my non-overlapping LC4K bed file seemed consistent with my LC_ProbesInfo file, (while Dr. Zhang's non-overlapping LC4K bed file was not, so I decided to continue with my non-overlapping LC4K bed file.
1631 ./bisPileupPair2specificity.pl LC_ProbesInfo_nonOverlappingBED.txt /media/TmpStore1/Noi_scratch/LC_20111017/mapping_hg18/H1-52C.fwd.pileup /media/TmpStore1/Noi_scratch/LC_20111017/mapping_hg18/H1-52C.rev.pileup > H1-52C.LC.specificity & 1632 ./bisPileupPair2specificity.pl LC_ProbesInfo_nonOverlappingBED.txt /media/TmpStore1/Noi_scratch/LC_20111017/mapping_hg18/H1-54C.fwd.pileup /media/TmpStore1/Noi_scratch/LC_20111017/mapping_hg18/H1-54C.rev.pileup > H1-54C.LC.specificity & 1633 ./bisPileupPair2specificity.pl LC_ProbesInfo_nonOverlappingBED.txt /media/TmpStore1/Noi_scratch/LC_20111017/mapping_hg18/PGP1F-54C.fwd.pileup /media/TmpStore1/Noi_scratch/LC_20111017/mapping_hg18/PGP1F-54C.rev.pileup > PGP1F-54C.LC.specificity & 1634 ./bisPileupPair2specificity.pl LC_ProbesInfo_nonOverlappingBED.txt /media/TmpStore1/Noi_scratch/LC_20111017/mapping_hg18/PGP1F-52C.fwd.pileup /media/TmpStore1/Noi_scratch/LC_20111017/mapping_hg18/PGP1F-52C.rev.pileup > PGP1F-52C.LC.specificity & 1639 ./bisPileupPair2specificity_DD.pl LC_ProbesInfo_nonOverlappingBED.txt /media/TmpStore1/Noi_scratch/LC_20111017/mapping_hg18/PGP1iPS-52C.fwd.pileup /media/TmpStore1/Noi_scratch/LC_20111017/mapping_hg18/PGP1iPS-52C.rev.pileup > PGP1iPS-52C.LC.specificity & 1640 ./bisPileupPair2specificity_DD.pl LC_ProbesInfo_nonOverlappingBED.txt /media/TmpStore1/Noi_scratch/LC_20111017/mapping_hg18/PGP1iPS-54C.fwd.pileup /media/TmpStore1/Noi_scratch/LC_20111017/mapping_hg18/PGP1iPS-54C.rev.pileup > PGP1iPS-54C.LC.specificity &
- H1-54C
Size of target = 257866 Total mappable bases = 168941334 Total target bases = 89608268 (0.530) Total target covered 1X = 185056 (0.718)
- H1-52C
Size of target = 257866 Total mappable bases = 189704106 Total target bases = 78915250 (0.416) Total target covered 1X = 187795 (0.728)
- PGP1iPS-54C
Size of target = 257866 Total mappable bases = 210319026 Total target bases = 131291469 (0.624) Total target covered 1X = 185174 (0.718)
- PGP1iPS-52C
Size of target = 257866 Total mappable bases = 257710476 Total target bases = 124280612 (0.482) Total target covered 1X = 188753 (0.732)
- PGP1F-54C
Size of target = 257866 Total mappable bases = 301811294 Total target bases = 159092545 (0.527) Total target covered 1X = 188938 (0.733)
- PGP1F-52C
Size of target = 257866 Total mappable bases = 301526024 Total target bases = 121979569 (0.405) Total target covered 1X = 193326 (0.750)
1645 ./bisPileupPair2specificity_DD.pl CpG330KProbesInfo_nonOverlappingBED-10292011.txt UCLA_GK0027-004.fwd.pileup UCLA_GK0027-004.rev.pileup > UCLA_GK0027-004.330K.specificity &
Size of target = 35827535 Total mappable bases = 1687630717 Total target bases = 1582327158 (0.938) Total target covered 1X = 25050732 (0.699)
1647 ./bisPileupPair2specificity_DD.pl CpG330KProbesInfo_nonOverlappingBED-10292011.txt /media/Ext9T/DD_Ext9T/KangZhang_PE110/11201001c.txt.fwd.pileup /media/Ext9T/DD_Ext9T/KangZhang_PE110/11201001c.txt.rev.pileup > KangZhang.11201001c.330K.specificity &
Size of target = 35827535 Total mappable bases = 2340766083 Total target bases = 2254282508 (0.963) Total target covered 1X = 24118127 (0.673)
- HL094 PGP1F
Size of target = 35827535 Total mappable bases = 1019718323 Total target bases = 934477017 (0.916) Total target covered 1X = 13888586 (0.388)
- HL094 H1
Size of target = 35827535 Total mappable bases = 1292246322 Total target bases = 1150312135 (0.890) Total target covered 1X = 8637207 (0.241)
- HL094 PGP1iPS
Size of target = 35827535 Total mappable bases = 1900359617 Total target bases = 1830971383 (0.963) Total target covered 1X = 25334157 (0.707)
Discussion[edit]
- 1 LC4K appears to be more sensitive (>70% target base covered at >=1X). But it is less specific.
- 2 In the recent whole blood capture with CpG330K has a high specificity (>90% mappable bases were on target), but still ~30% of target is missing. Are the similar sites missing between different samples?
- 3 Maybe in general, every genome will be missing or have very difficult to capture segments at a rate of ~30% (of our probes set.) This could be because of nearby SNPs or unconserved part of the genome. It'll be interesting to see what 30% are keep getting missed.