Brandon:LabNotes/Project1/2012-7-27

From ZhangLabWiki
Jump to navigation Jump to search

Peak calling and data analysis, continued from IVT method accessibilty comp[edit]

  • 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[edit]

  • issues found in bam files, bed files, using macs14, peak calling
  • correct input data for MACS14 is unique BAM files OR BED files created by using the bedtools function bamToBed.
  • can use genomeCoverageBed to make files to view how much reads there are on UCSC genome browser since smaller files are needed
  • can use MACS generated .bed file to see significant region in a single bar, but cannot see approximate amount of reads per region.


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. Ran MACS peak calling 3 times on each type of data to ensure reproducible.
    • 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


  • more peaks called when using bam files. Less BP called when using BAM files. Not much difference seen in control samples for peaks called and especially for BP called.

File:Peaks statistics.png


  • number of base pairs called from the different files. MACS calls less sgifnigicant basepairs from the bam files as opposed to the bed extracted files w/ awk.
' bed_extracted bam bed_coverage '
bp_called2 bp_run2 bp_called2 % difference
10000 185,726 145,647 1,313,689 0.784203612
1000 1,403,020 510,511 5,343,351 0.363865804
100 63,507 56,258 271,512 0.885855103
6ng 345,760 317,124 3,182,917 0.917179547
600pg 31,723 34,255 193,392 1.079815906
60pg 6,033 6,147 7,100 1.018896072

File:Peaks statistics2.png


Correct way to call peaks is with BAM files or BED from created from using bamToBed

  • generated a BED file using bamToBed and compared it to peaks and basepairs called when using a unique BAM file (after sorting and clonal read removal)
' bamToBed bam bamToBed bam
peaks_run2 peaks_run2 bp_run2 bp_run2
10000 656 656 145,646 145,647
1000 3,343 3,343 510,510 510,511
100 360 360 56,258 56,258
6ng 564 564 317,124 317,124
600pg 20 20 34,255 34,255
60pg 24 24 6,147 6,147




side by side data comparison of nextera accessibility assay versus IVT assay[edit]

  • issue of clonal reads has been addressed with the IVT method. However now have to improve efficiency of the assay by reducing the number of failed reads.


first two columns are IVT, second two columns are nextera. different cell numbers used but still can draw conclusions

sample clonal read% sample clonal read%
10000 cells 1.982% gentle-lysis-500-cell 89.22%
1000 cells 4.176% gentle-lysis-200-cell 89.34%
100 cells 6.158% gentle-lysis-100-cell 89.49%
6 ng pure DNA 2.131% lyse+protease-500-cell 5.58%
600 pure pg 1.159% lyse+protease-200-cell 5.38%
60 pure pg 0.468% lyse+protease-100-cell 7.41%
sample unique reads sample unique reads
10000 cells 1,116,444 gentle-lysis-500-cell 862,574
1000 cells 1,551,797 gentle-lysis-200-cell 290,269
100 cells 332,102 gentle-lysis-100-cell 200,865
6 ng pure DNA 2,924,988 lyse+protease-500-cell 3,508,328
600 pure pg 392,021 lyse+protease-200-cell 2,912,333
60 pure pg 21,689 lyse+protease-100-cell 1,994,719
sample fail map % sample fail map %
10000 cells 65.52% gentle-lysis-500-cell 22.07%
1000 cells 73.89% gentle-lysis-200-cell 25.03%
100 cells 91.80% gentle-lysis-100-cell 36.30%
6 ng pure DNA 58.66% lyse+protease-500-cell 9.98%
600 pure pg 86.91% lyse+protease-200-cell 11.90%
60 pure pg 98.99% lyse+protease-100-cell 12.54%


  • comparions of peaks called in different sample types.
  • from file _peaks_calling_data_analysis_etc.xlsx in data analysis folder


  • less peaks called in pure/protease treated DNA, which is good since should be more random
  • similar amount of peaks called with IVT and nextera data. but nextera tagmentation data seemed to be of lower quality than the IVT data.
  • UCSC Genome browser data (run though my own analysis, same as other data)


UCSC bam bam
UCSC_peaks UCSC_bp_cvd
rep1 91,965 9,470,220
rep3 107,754 19,738,155
rep4 107,754 14,323,016
rep5 122,159 17,820,762
average 107,408 15,338,038



  • IVT assay


IVT bam bam
IVT_peaks IVT_bp_cvd
10000 656 145,647
1000 3,343 510,511
100 360 56,258
6ng 564 317,124
600pg 20 34,255
60pg 24 6,147



  • nextera tagmentation assay


nextera bam bam
nxta_peaks nxta_bp_cvd
500-cell-acc 5,698 403,529
200-cell-acc 258 131,013
100-cell-acc 794 259,172
500-cell-pro 1,705 164,349
200-cell-pro 865 89,491
100-cell-pro 574 63,078


  • similar amount of peaks called with IVT and nextera data. but nextera tagmentation data seemed to be of lower quality than the IVT data.


  • looked at how well data mapped to the UCSC genome browser GM12878 dnase I hypersensitivity tract for all the samples. Used coverageBed function.
  • calculated overlap by overlaying UCSC data onto IVT data, where any sequence called in UCSC but not called in the IVT data is ignored. If IVT data was overlaid onto UCSC data, there would only be 2% or less overlap since so few bases are called in IVT data versus UCSC data.


  • important information on how coverageBed analysis was performed: samples with larger number number of basepairs called in peaks were always overlaid onto samples where less basepairs were called in peaks. Exception is when the analysis was done both ways with a sample regardless of how many basepairs were called in the peaks. For example the 1000cellIVT (510,511 bp called) sample is overlaid onto the 10000cellIVT (145,647) sample. That way any bp called in the 510,511 bp and are not called in the 145,647 sample are not calculated in the analysis.
    • same for UCSC genome browser data (19.7 million bp called) versus any IVT or nextera tagmentation sample called (most was 510,511 bp)
  • raw data and calculations after performing coverageBed
coverageBed overlap_bp total_bp % covearge '
UCSCRep3_onto_IVT1000cell 35,694,646 93,107,820 38.337% before_macs_peak_calling
UCSCrep3_onto_IVT10000cells 75,319 145,647 51.713% macs_called_peaks_bp_basis
UCSCrep3_onto_IVT1000cells 364,391 510,511 71.378% macs_called_peaks_bp_basis
UCSCrep3_onto_IVT100cells 30,247 56,258 53.765% macs_called_peaks_bp_basis
IVT1000cell_onto_UCSCrep3 364,391 26,352,152 1.383% macs_called_peaks_bp_basis
UCSCRep3_onto_nextera500cellacc 12,397,567 31,052,664 39.924% before_macs_peak_calling
UCSCrep3_onto_nextera500cellacc 63,106 403,529 15.639% macs_called_peaks_bp_basis
UCSCrep3_onto_200cellacc 11275 131013 8.606% macs_called_peaks_bp_basis
UCSCrep3_onto_nextera100cellacc 27310 259172 10.537% macs_called_peaks_bp_basis
UCSCrep3_onto_500cellpro 23,750 164,349 14.451% macs_called_peaks_bp_basis
UCSCrep3_onto_IVT6ngpure 16,655 317,124 5.252% macs_called_peaks_bp_basis
500cellacc_onto_500cellpro 24,590 164,349 14.962% macs_called_peaks_bp_basis
500cellPro_onto_500cellacc 24,590 403,529 6.094% macs_called_peaks_bp_basis
200cellacc_onto_500cellPro 11,034 164,349 6.714% macs_called_peaks_bp_basis
500cellpro_onto_200cellacc 11,034 131,013 8.422% macs_called_peaks_bp_basis
IVT1000cell_onto_IVT6ng 14,567 317,124 4.593% macs_called_peaks_bp_basis
IVT6ng_onto_IVT1000cell 14,567 510,511 2.853% macs_called_peaks_bp_basis
IVT10000cells_onto_IVT100cell 12,863 56,258 22.864% macs_called_peaks_bp_basis
IVT1000cell_onto_IVT10000cell 35,599 145,647 24.442% macs_called_peaks_bp_basis
IVT1000cell_onto_IVT100cell 17,820 56,258 31.675% macs_called_peaks_bp_basis
500cellacc_onoto_200cellacc 16,582 131,013 12.657% macs_called_peaks_bp_basis
500cellacc_onto_100cellacc 16,816 259,172 6.488% macs_called_peaks_bp_basis


  • overlap of UCSC peaks with called peaks. Can see nextera tagmentation data did not have strong overlap, but IVT data did. Though nextera tagmentation data was generated with GM20431 cells, but is also a lymphocyte line so should have similar peaks. hard to conclude if difference is only because they are different cell lines.

File:Picture1uyghgf.png


  • Overlap of UCSC peaks onto controls peaks (pure DNA or protease treated DNA before performing tagmentation). with nextera tagmentation data, 14% overlap is seen, which is the same as the accessibility samples. Thus either same peaks were called or accessibility assay was not as effective.
  • UCSC data onto control for IVT assay shows nearly 10% less peak overlap, thus more random. Though in this case, pure DNA was a IPS cell line, which could account for the difference. However UCSC data is difference lymphocyte line (GM12878) than nextera tagmentation data (GM20431)

File:Picture1.png


  • looking at overlap between accessibility samples and controls. Both nextera tagmentation data (columns 1-4) and IVT data (columns 5-6) compared.
  • IVT data still has less overlap between with controls than does the nextera tagmentation data. Also has 14% overlap again, so possibly the same universally accessible peaks were called in either case.

File:Picture3.png


  • next looked at reproducability between samples.
  • again IVT samples have much higher reproducibilty between the samples than does the nextera tagmentation data. This lends more evidence to the IVT assay having higher accuracy and reproducablilty overall, when compared to the nextera tagmentation assay.

File:Picture4.png

re-analysis of unmappable sequenning data[edit]

  • looked at reads where there is only a small insert or no insert at all.
  • seems to be added when random nonamer is added since sequence that is in random nonamer is present in ~400/470 of the reads. Thus is isnt added only by the illuminia PCR primer thats for the 3' end.
  • analyszed the file 2012-07-17 unmappble index 51 reads pure dna.docx in the sequencing folder for this sequening run


  • also checked if could possibly be a concatamer. doesn't seem like it.
    • random nonamers could be binding to each other but then shouldn't be singing complete adatptor/illuminia adaptor ends, which is what is seen. Then they could have a random nonamer sequence in the middle of of the adaptor end, and the start of another adaptor. but that is not seen.
    • also some 3' adatptors are cut off and thus the whole illuminia adaptor and piece of 3' adaptor are sequenced. The end of the illuminia adatpors is reached and a series of A's are sequenced. Thus the true end has been reached, otherwise if it was concatenated then the start of another 3' adaptor/illuminia adaptor would be seen.
  • Most sequences for nonamer primed RT synthesis are either adajacent to the 5' sequence
  • For RNase III fragmented RNA the same is seen. Transposons containing the 5' read sequence are still left over. They are made into 2x stranded DNA with taq polymerase extension, then during IVT can generate small RNA fragments consisting only of the read sequence. Next, this fragment is polyadenylated with polyA polymerase, and thus a fragment containing only the 5' end, and polyA is created. subsequent amplification results in a insertless fragment.


  • To address these issues, Rnase III method/PAP polyA tail can be used instead of random nonamer. Or Mg++ fragmentation and end repair, though requires more steps.
  • Also cleaning of RNA with Zymo clean and concentrator kit and performing DNase I cleaning before further steps. Make sure to use Zymo protocol removing fragments less than 200 bp.
  • should results in alot less unmappable reads.


conclusions[edit]

  • seems that the IVT assay has higher accurarcy and reproduceability than the nextera tagmentation assay after doing the coverageBed analysis. There is still much optimization that can be done with the IVT assay too, while not much can be done with the nextera assay without major modifications.
    • similar number of peaks were called with IVT assay, but were of higher quality than the nextera tagmentation assay.


  • peak calling will be done using the right files now.
  • optimizations include:
    • Dnase I digestion after IVT
    • running Zymo RNA purification so that fragments < 200 bp will be removed. (before ran so anything under 17 would be removed.
    • using RNase III fragmentation/PolyA addition opposed to RT with random nonamers