Dinh 2011/NOTES/2011-10-29

From ZhangLabWiki
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
./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.
./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. 1 LC4K appears to be more sensitive (>70% target base covered at >=1X). But it is less specific.
  1. 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?
  1. 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.