Brandon:LabNotes/Project1/2012-10-19
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.
- made powerpoint slides of results, and made graphic of the IVT method detailing the method.
- need to fix up UCSC genome browser overlay pictures so they are more readable and and text is easier to see. also so font is the same in each of the pix.
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");
- bed to wig conversion script used (bed2wig.pl)
- script persent here: bed file of reads to wiggle file of numbers of features at each position
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.
- 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.
- 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