Brandon:LabNotes/Project1/2012-7-27
Jump to navigation
Jump to search
Peak calling and data analysis, continued from IVT method accessibilty comp
- seems to be working. comparing to GM12878 data to see effectiveness
- macs seems to be calling more significant peaks when analyzing from extracted bed files, but less total bps in covered areas.
- analyzed previous data wrong since data in bed files did not overlap.
- thus if 2 reads overlapped, it would be 3 counts in the bed file.
- if 3 did it would count as 5 reads
- if 6 did, it would be read as ~10 reads.
bam, bed issues, peak calling discrepancies
- issues found in bam files, bed files, using macs14, peak callind
peak calling issues BED files
- when using bedtools to call genome coverage with the coverage.bed function, it results in formation of bed entries that to not overlap. Thus if 3 reads partially overlap, 5 entries will be made based on the points where reads start/begin
genomeCoverageBed -bg -ibam $uniqueBam -g /home/bsos/software/BEDTools-Version-2.12.0/genomes/human.hg19.genome > $bed
- This results in artificially inflated read counts for for overlapping reads, which is why many more peak calls are made.
- Though this is useful for viewing on UCSC genome browser which requires smaller files as inputs.
- all analyses are being fixed by not using BEDs generated from using genomeCoverageBed. (use SAM files, or extracted BED files from SAM files)
- 6 reads reported from BAM file (IVT_extracted_BED tract) are reported as 10 reads in the genomeCoverageBed generated file (IVT_1000_cell tract).
File:BAMorExtractedBed versus coverageBED.png
- 3 reads reported from BAM file (IVT_extracted_BED tract) are reported as 5 reads in the genomeCoverageBed generated file (IVT_1000_cell tract).
File:BAMorExtractedBed versus coverageBED2.png
bam/extracted bed file calling discrepancies
- macs14 seems to be better at calling peaks directly from the bam file. THus BED files should probabaly be only reserved for viewing reads etc on the UCSC genome browswer. As can be seen, more peaks are called in the accessibility assay on cells versus pure DNA, which is what we should be seeing.
- first column is extracted BED file from BAM.
samtools view s_5_1_Indx49.unique.bam | awk '{print $3"\t"$4"\t"$4+60}'| sort -u > s_5_1_Indx49.unique.bed
- second columns is unique bam file
- 3rd column is bed file from using genomeCoverageBed
' | bed_extracted | bam | bed_coverage | ' |
peaks_run2 | peaks_run2 | peaks2 | % difference | |
10000 | 348 | 656 | 3143 | 1.885057471 |
1000 | 2,549 | 3,343 | 11247 | 1.311494704 |
100 | 131 | 360 | 651 | 2.748091603 |
6ng | 642 | 564 | 6529 | 0.878504673 |
600pg | 60 | 20 | 447 | 0.333333333 |
60pg | 24 | 24 | 25 | 1 |
all | most | all | ||
- check UCSC data and see how that compares. make graphs?
- ran replicates to ensure peak calling was similar between runs and was reproducible.