Brandon:LabNotes/Project1/2012-10-19

From ZhangLabWiki
Jump to navigation Jump to search

Continuation/Update to Accessability DNA sequence analysis from IVT Amplification round 2 after more analysis[edit]

  • combined new sequencing data run to draw further conclusions
  • downsampled samples so similar number of reads were analyzed
  • did comparison and statistics analysis of data between IVT data and nextera data
  • created scripts to extract specific regions from bam files so can easily view accessability and compare samples in UCSC genome browser. uploading whole genome accessabilty data to the UCSC genome browser is not easily feasible for what needs to be accomplished.



scripts used for extracting a genomic region and converted from bam to wig to visualize read densities[edit]

  • so can be visualized on the UCSC genome browser under custom tracts.


  • need a bam file of analyzed read results.
  • Used unique.bam files from my pipeline. (are sorted and have clonal reads removed)


  • outline: bam file -> index bam file -> extract region -> convert to bed -> convert to wig -> add header -> concatenate wig files
  • wig files are concatenated so only need to be uploaded to the UCSC genome browser in one file.


  • shell script, GO.runextractBamToWig.sh
  • only need to modify argument 2, chromosome number and genomic location, which is the region you want to extract. Then run the shell script to extract the region you want.
for f in *.bam
do
	perl /media/Ext12T/BS_Scratch2/121004_HL131/UCSC_peak_overlap_bam_extraction_to_wig/extractBamToWig.pl $f chr19:27,711,974-27,752,373
done;


cat *ded.wig > all_combined.wig
echo -e "\n *ded.wig > all_combined.wig \n"


  • perl script used: extractBamToWig.pl
#!/usr/bin/perl -w

#input:   perl testing.pl dc_Indx50.unique.bam chr21:9822000-9832000

use strict;

my $in_file = $ARGV[0];
my $chr_location = $ARGV[1];

my $split = $in_file;
$split =~ s/.bam//;  #substitution

my $extracted_bam = $split . "_" . $chr_location . ".bam";
my $bed_file = $split . "_" . $chr_location . ".bed";
my $wig_file = $split . "_" . $chr_location . ".wig";
my $final_wig = $split . "_" . $chr_location . "hadded" . ".wig";
my $index_file = $in_file . ".bai";


my $track_header = qq{browser position $chr_location\ntrack type=wiggle_0 name="$split" description="$split"  visibility=full autoScale=off viewLimits=0.0:25.0 color=50,150,255 yLineMark=11.76 yLineOnOff=on priority=10\n };

open (MYFILE, '>>header.txt');
print MYFILE "$track_header";
close (MYFILE);

#print "\n$track_header\n";
#print "\n$in_file\n";
#print "\n$split\n";
#print "\n$chr_location\n";

my $cmd = "samtools index $in_file";
print "\n$cmd\n";
system($cmd);


$cmd = "samtools view -bh $in_file $chr_location > $extracted_bam";
print "\n$cmd\n";
system($cmd);


$cmd = "bamToBed -i $extracted_bam > $bed_file";
print "\n$cmd\n";
system($cmd);


$cmd = "perl bed2wig.pl $bed_file > $wig_file";
print "\n$cmd\n";
system($cmd);


$cmd = "cat header.txt $wig_file > $final_wig";
print "\n$cmd\n\n";
system($cmd);

unlink("header.txt");
unlink("$wig_file");
unlink("$bed_file");
unlink("$extracted_bam");
unlink("$index_file");


use strict;
use warnings;

my $file = $ARGV[0];
my %hash = ();

open(FILE, "<$file") || die "cannot open $file\n";
while(<FILE>){
    chomp;
    my ($chr, $start, $end, $id, $score, $strand) = split(/\s+/,$_);
    for(my $i=$start; $i<=$end; $i++){
        $hash{$chr}{$i}++;
    }
}
close(FILE);

foreach my $chr (keys %hash){
    print "variableStep chrom=$chr\n";
    foreach my $i (sort {$a<=>$b} keys %{$hash{$chr}}){
        print "$i\t$hash{$chr}{$i}\n";
    }
}


sequencing results[edit]

  • statistics and results from combining all the sequencing runs
  • downampled samples were used in all the following comparison analyses for comparing the IVT method to the Nextera method.
  • data combined from HL130, HL131, DD04 run
  • odd that there were more clonal reads in the IVT samples this round, not sure why. but 1000 cell sample worked pretty well. still called many more peaks in IVT samples


' sample total reads one alignment one alignment % failed to map fail map % supressed supressed % clonal reads clonal read% unique reads paired peaks peaks
1000 cells IVT s_7_1_Indx49 11,921,212 8,277,293 69.43% 2,131,831 17.88% 1,512,088 12.68% 308,039 3.72% 7,969,254 654 24,848
500 cells IVT s_7_1_Indx50 12,846,634 8,135,684 63.33% 2,544,528 19.81% 2,166,422 16.86% 4,213,575 51.79% 3,922,109 576 9,394
100 cells IVT s_7_1_Indx51 11,696,728 4,083,894 34.91% 2,626,316 22.45% 4,986,518 42.63% 2,213,554 54.20% 1,870,340 373 1,989
6 ng purified DNA IVT s_7_1_Indx52 20,394,388 13,419,300 65.80% 3,552,094 17.42% 3,422,994 16.78% 6,927,156 51.62% 6,492,144 102 3,112
600 pg purified DNA IVT s_7_1_Indx53 2,761,684 343,038 12.42% 702,766 25.45% 1,715,880 62.13% 198,583 57.89% 144,455 - 476
60 pg purified DNA IVT s_7_1_Indx54 603,428 27,006 4.48% 170,258 28.22% 406,164 67.31% 18,584 68.81% 8,422 - 315
1000 cells Nxtra s_7_1_Indx55 8,585,321 5,371,720 62.57% 1,989,397 23.17% 1,224,204 14.26% 2,541,231 47.31% 2,830,489 2,194 3,279
500 cells Nxtra s_7_1_Indx56 28,365,285 17,189,378 60.60% 7,478,535 26.37% 3,697,372 13.03% 1,490,960 8.67% 15,698,418 541 12,793
100 cells Nxtra s_7_1_Indx57 42,902,381 26,430,162 61.61% 11,099,059 25.87% 5,373,160 12.52% 6,319,086 23.91% 20,111,076 379 13,631
6 ng purified DNA Nxtra s_7_1_Indx58 23,668,232 13,653,302 57.69% 6,268,204 26.48% 3,746,726 15.83% 538,700 3.95% 13,114,602 168 4,455
600 pg purified DNA Nxtra s_7_1_Indx59 3,908,464 1,831,198 46.85% 1,343,607 34.38% 733,659 18.77% 939,230 51.29% 891,968 34,477 7,566
60 pg purified DNA Nxtra s_7_1_Indx60 1,967,066 893,437 45.42% 790,544 40.19% 283,085 14.39% 554,383 62.05% 339,054 8,974 988



  • downsample results
  • samples 49, 50, 52, 56, 57 58 downsampled to 100 cell IVt total reads amount which was the lowest one in the group above 10 million reads


' sample total reads one alignment one alignment % failed to map fail map % supressed supressed % clonal reads clonal read% unique reads paired peaks peaks
1000 cells IVT s_7_1_Indx49 11,686,667 8,114,449 69.43% 2,089,796 17.88% 1,482,422 12.68% 297,361 3.66% 7,817,088 591 24,674
500 cells IVT s_7_1_Indx50 11,677,578 7,395,709 63.33% 2,312,854 19.81% 1,969,015 16.86% 3,505,559 47.40% 3,890,150 543 9,583
100 cells IVT s_7_1_Indx51 11,696,728 4,083,894 34.91% 2,626,316 22.45% 4,986,518 42.63% 2,213,554 54.20% 1,870,340 373 1,989
6 ng purified DNA IVT s_7_1_Indx52 11,723,300 7,712,924 65.79% 2,042,896 17.43% 1,967,480 16.78% 2,369,254 30.72% 5,343,670 128 7,047
600 pg purified DNA IVT s_7_1_Indx53 2,761,684 343,038 12.42% 702,766 25.45% 1,715,880 62.13% 198,583 57.89% 144,455 - 476
60 pg purified DNA IVT s_7_1_Indx54 603,428 27,006 4.48% 170,258 28.22% 406,164 67.31% 18,584 68.81% 8,422 - 315
1000 cells Nxtra s_7_1_Indx55 8,585,321 5,371,720 62.57% 1,989,397 23.17% 1,224,204 14.26% 2,541,231 47.31% 2,830,489 2,194 3,279
500 cells Nxtra s_7_1_Indx56 11,816,738 7,161,353 60.60% 3,116,151 26.37% 1,539,234 13.03% 330,186 4.61% 6,831,167 886 4,722
100 cells Nxtra s_7_1_Indx57 11,687,750 7,201,722 61.62% 3,023,309 25.87% 1,462,719 12.51% 627,951 8.72% 6,573,771 674 4,115
6 ng purified DNA Nxtra s_7_1_Indx58 11,384,161 6,568,137 57.70% 3,015,484 26.49% 1,800,540 15.82% 139,063 2.12% 6,429,074 195 2,324
600 pg purified DNA Nxtra s_7_1_Indx59 3,908,464 1,831,198 46.85% 1,343,607 34.38% 733,659 18.77% 939,230 51.29% 891,968 34,477 7,566
60 pg purified DNA Nxtra s_7_1_Indx60 1,967,066 893,437 45.42% 790,544 40.19% 283,085 14.39% 554,383 62.05% 339,054 8,974 988


macs mfold testing[edit]

  • tested different mfold values in macs14, since they recommended to try different settings if the paired peaks values were less than 1000. mfold sets the stringency for cut offs and inclusion/exclusion of peaks.


  • Tried different settings, result in /media/Ext12T/BS_Scratch2/121004_HL131/macs_mfold_testing on genome minor


  • found that basically there was no change on the results when mfold values were changed so that cutoffs were less stringent. more paired peaks were called but the same number of normal peaks/basepairs were called.


  • when mfold values were more stringent, as expected less peaks were called, and less bp were covered.


  • Thus will not change mfold values when running macs14 for peak calling





statistics showing IVT method working and better than nextera[edit]

  • comparative metrics used:
  • Overlap of accessible regions called in UCSC genome browser GM12878 DNase hypersensitivity data, based on number of reads needed.


  • Visualization of read densities with new methods compared to already called accessible regions on the UCSC genome browser


  • Number of peaks found with analysis by MACS14


  • Clonal reads, Unique reads, informative data obtained from read amounts and coverage.


  • Analysis of pure DNA controls



  • Overlap of peaks called when data processed by MACS14.
  • Used Bedtools, coverageBed function for analysis. Ran UCSC genome browser GM12878 data though my pipeline. Called peaks on UCSC genome browser with MACS14. Overlaid UCSC genome browser data onto peaks called in IVT data.
  • Nextera has affinity for same regions, in purified DNA versus untreated DNA
  • Bias to specific regions for unknown reasons
  • Technically percent overlap for pure DNA controls should be 1.35% since accessible BP in UCSC genome browser is 40,587,687 out of 3 billion BP.
  • Adds noise to the analysis or could be over representing accessible regions it has an affinity for.

File:Picture1.png


  • Performed with coverageBed, overlap onto UCSC genome browser
  • Need 2.8X more unique reads in nextera to obtain the same number of accessible BP covered with the IVT assay.
  • Need 2.18X more total reads in nextera to obtain the same number of accessible BP covered with the IVT assay.
  • First is UCSC data overlaid onto IVT data, then NXTA data for the second.
UCSC onto all combined for combined IVT or NXTA ' ' ' ' '
not downsampled
% overlap UCSC unique reads ttl total reads before analysis
UCSC_all_onto_IVT_all.txt 2,583,808 4,279,514 0.603762016 13,761,703 36,464,574
UCSC_all_onto_NXTA_all.txt 3,224,000 4,393,788 0.733763213 38,639,983 79,852,987
2.807790794 2.189878511



  • Performed with coverageBed, overlap onto UCSC genome browser
  • Downsampled data so nextera is comparable to IVT
  • Unique read count is similar, as is total reads.
  • IVT method predicts 3.36X more accessible regions than the nextera method with the same amount of total reads and unique reads.
UCSC onto all combined for combined IVT or NXTA ' ' ' ' ' '
downsampled unique reads total reads before analysis
UCSC_all_onto_IVT_combined.txt 2,662,887 4,447,491 0.598739154 13,577,578 35,060,973
UCSC_all_onto_NXTA_combined.txt 871,396 1,324,236 0.658036785 16,235,427 32,089,809
3.358533524 3X more efficient?



peaks tables

  • IVT samples only
' total reads unique reads peaks called
1000 cells IVT 11,921,212 7,969,254 24,848
500 cells IVT 12,846,634 3,922,109 9,394
100 cells IVT 11,696,728 1,870,340 1,989
6 ng purified DNA IVT 20,394,388 6,492,144 3,112
600 pg purified DNA IVT 2,761,684 144,455 476
60 pg purified DNA IVT 603,428 8,422 315



  • in comparison with nextera
sample peaks called
1000 cells IVT 24,674
1000 cells NXTA 3,279
500 cells IVT 9,583
500 cells NXTA 4,722
100 cells IVT 1,989
100 cells NXTA 4,115


  • Overall will need less sequencing depth with IVT method to obtain the same data
  • Total BP called from UCSC genome browser data: 40,587,687 and Total peaks: 181,595
  • Collected data and plotted peaks called versus number of reads from analyses
  • Extrapolation suggests 125 million unique reads needed for IVT and 492 million unique reads needed for comparable coverage represented by the UCSC data.

File:Picture2.png


  • Overall data generated from IVT samples is better able to call peaks. Data from nextera seems to be dependent on the number of reads, regardless of cellular concentration, not sure why.


  • Also seems like nextera tagmentation might have a bias for the same regions, since pure DNA samples do less of a random distribution than IVT samples.


Comparison data on UCSC genome browser[edit]

  • check powerpoint slides up in top, at this point those need to be redone to look presentable.


blah[edit]

  • bed
  • look at unmappables
  • plan next experiemntsd for RNA-seq and accessibility