Brandon:LabNotes/Project1/2012-7-27: Difference between revisions
Jump to navigation
Jump to search
>Bsos |
>Bsos |
||
Line 223: | Line 223: | ||
*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. | *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 | *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. |
Revision as of 22:19, 2 August 2012
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 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.
- 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 |
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
- 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.
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% |
- finish updating wiki on reads eetc
- check UCSC data and see how that compares. make graphs?
- ran replicates to ensure peak calling was similar between runs and was reproducible.
differences in peaking calling nbtn extracted bed/bam
comparision between GM12878 datas
re-analysis of unmappable sequenning data
- 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.