Hi Alex, your pipeline is almost perfect for my need, except that the resulted bed doesn't contain the strand information, which can be easily done with additional parameter to your awk code. Do you think it's possible to add counts of the duplicate reads in resulted bed file?
annotation of duplicate reads
Hi there,
I have a lot of duplicate reads in my RNA-seq, which is great from my view. My question is is there any way to get the coordinates of these duplicate reads (> 10 reads) from a BAM/SAM file? It will be perfect to be able to annotate these regions with overlapping or flanking genes and output the results in a bed format. Any suggestions are appreciated, but solutions using R packages are more preferred.
Here are some artificial example datasets (might be deleted after one year) to be tested.
2 answers
A lot of this can be done easily and quickly with the BEDOPS toolkit and standard UNIX tools:
To convert BAM to BED:
$ bam2bed < reads.bam > reads.bed
To get a list of duplicate reads that meet your 10-read threshold, map the reads to themselves with bedmap, using its --count operator with the --exact overlap threshold:
$ bedmap --count --exact --echo --delim '\t' reads.bed \
| awk '$1>=10' \
| cut -f2- \
> duplicate_reads.bed
To get a list of the genomic regions of the duplicate reads:
$ awk 'BEGIN {pChr=""; pStart=-1; pStop=-1; } \
{ \
cChr=$1; cStart=$2; cStop=$3; \
if ((pChr!=cChr) || (pStart!=cStart) || (pStop!=cStop)) { \
print $1"\t"$2"\t"$3; \
pChr=cChr; pStart=cStart; pStop=cStop; \
} \
}' duplicate_reads.bed > duplicate_regions.bed
To annotate these, we can use gtf2bed to convert GENCODE v20 protein-coding elements to a BED file, bedops to pad duplicate regions (say, a 1000-base window), and bedmap to map the GENCODE gene annotations to the duplicate regions.
To get some annotations of interest into BED format, for example:
$ wget -qO- ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_20/gencode.v20.annotation.gtf.gz \
| gunzip -c - \
| gtf2bed - \
| grep -i protein_coding - \
> gencode.v20.pcRegions.bed
To pad regions and map these GENCODE v20 protein-coding annotations to those padded regions:
$ bedops --everything --range 1000 duplicate_regions.bed \
| bedmap --echo --echo-map-id-uniq - gencode.v20.pcRegions.bed \
> annotated_and_padded_duplicate_regions.bed
• 0 views
•
link
updated
by
Ram
•
written
by
Alex Reynolds
•
written
by
Alex Reynolds
The reads contain directionality, so you could modify the awk processing step:
$ awk 'BEGIN {pChr=""; pStart=-1; pStop=-1; pStrand="."; } \
{ \
cChr=$1; cStart=$2; cStop=$3; cStrand=$6; \
if ((pChr!=cChr) || (pStart!=cStart) || (pStop!=cStop)) { \
print cChr"\t"cStart"\t"cStop"\tid-"NR"\t.\t"cStrand; \
pChr=cChr; pStart=cStart; pStop=cStop; pStrand=cStrand; \
} \
}' duplicate_reads.bed > duplicate_regions.bed
This gives you a six-column, stranded BED file specifying duplicate regions.
To get counts, initialize a counter variable to zero in the BEGIN block. When you find a BED line with different records, you could modify this procedure to print the previous BED element's chromosome, start and stop value, along with the value of the counter that you keep iterating. One you have printed some element, you reset counter to zero. In other words, instead of printing cChr etc., you would print pChr, etc. and counter as the score value.
Finally, you would add an END block (not shown) that prints the state at the very last line of the script, if the last line has the same coordinates as the previous-to-last-line. That makes sure you print everything.
It's a slightly more complex awk script, as you have to keep track of a little more state, but this description should get you started.
To do strand-specific operations, you'd split the regions by directionality:
$ awk '$6=="+"' duplicate_regions.bed > duplicate_regions.sense.bed
$ awk '$6=="-"' duplicate_regions.bed > duplicate_regions.antisense.bed
You'd also need to split the GENCODE annotations by directionality:
$ awk '$6=="+"' gencode.v20.pcRegions.bed > gencode.v20.pcRegions.sense.bed
$ awk '$6=="-"' gencode.v20.pcRegions.bed > gencode.v20.pcRegions.antisense.bed
In the last step where you map GENCODE annotations to padded regions, you would do two bedmap operations, one on each of the directions for the regions reference file and annotations map file.
If you want the final answer in one file:
$ bedops --everything annotated_regions.sense.bed annotated_regions.antisense.bed > all_annotated_regions.bed
One complication I could see is in padding, but bedops --range with one value pads symmetrically, so that is taken care of.
• 2 views
•
link
updated
by
Ram
•
written
by
Alex Reynolds
•
written
by
Alex Reynolds
Change:
if ((pChr!=cChr) || ((pStart!=cStart) && (pStop!=cStop)))
To:
if ((pChr!=cChr) || (pStart!=cStart) || (pStop!=cStop))
(I'll edit my answers to fix this bug.)
• 2 views
•
link
written
by
Alex Reynolds
Duh. I just realized that you can reuse your duplicate reads, to map them back to duplicate regions and get your counts:
$ bedmap --echo --count --exact --delim '\t' duplicate_regions.bed duplicate_reads.bed \
| awk '{print $1"\t"$2"\t"$3"\tid-"NR"\t"$4}' \
> duplicate_regions_with_counts.bed
Or:
$ bedmap --echo --echo-map-id-uniq --count --exact --delim '\t' duplicate_regions.bed duplicate_reads.bed > duplicate_regions_with_ids_and_counts.bed
• 2 views
•
link
written
by
Alex Reynolds
Hi Alex, your pipeline worked perfectly on my example data set, but it seemed not efficient when dealing with real data. Even for a 100MB BAM input file, it took more than 24 hours to the 2nd step (get the duplicate_reads.bed) with my 4GB RAM Mac. Now it is running on a server (32 cores, 64 GB RAM) more than 1 hours, and I still don't see any hope.
I think I still need a efficient way, because there are more than 10 BAM files need to be analysed.
The conversion scripts are indeed slow and will be the bottleneck. I am working on a Python-free replacement of conversion tools for release in v2.5 (hopefully). In the meantime, you might test replacing bam2bed with bedtools bamToBed and see if it converts all the columns you need. It runs faster than v2.4.2 bam2bed. The rest of the pipeline should be very efficient.
• 0 views
•
link
written
by
Alex Reynolds
Log in to answer this question.
More posts like this
-
Customisable DNA Read Simulation for Deleted Reads
written by chris.fearn •
-
Annotating peaks with genes
written by V
-
annotation of chipseq peaks
written by rob.costa1234
-
How to extract best-hit reads from merged sam file with multiple hits per read?
written by O.rka
-
pancreatic cancer datasets
written by boymin2020 •
-
Extracting Codons from a SNP alignment
written by tpaisie
-
Join mapped, overlapping, paired-end reads
written by Mr. Dave •
-
Cluster the intervals but keep strand and return the average score
written by ifreecell
-
Count 5'End Mapped To A Specific Genomic Position
written by ifreecell
-
Duplicates In The Context Of Deep Sequencing
written by KCC
What have you tried? It'd be possible to do this in any language with convenient access to BAM files (C, python, R, java, etc.).
Could you just give some clue?
The general structure would be to read in a coordinate sorted BAM file and then create a buffer in a language of your choosing to hold one alignment. You'd then iterate through the alignments checking if a given alignment matches that in the buffer. If it does, increment a counter. If the alignment doesn't match that in the buffer, then look at the counter and see if there are enough duplicates to note the location. This would then be repeated. This is easy enough with single-end reads, but can be a bit more complicated with paired-end reads, depending on how you want to define duplicates.
If you are only interested in retrieving list of genes which have lots of duplicate reads then just compare the read count files generated using bam file before and after removal or flagging of duplicate reads. Read count softwares like HTseq will consider overlapping reads towards the counts. You can then compare these two count files and extract genes that show reduction of more than 10 counts between these two files. Retrieving the intergenic regions that show duplication is not that straightforward but can be accomplished using bedtools.