Greetings,
I have an RNA-seq dataset (HiSeq 3000, poly-A selected, 100 bp PE reads) with very high duplication, and I want to create coverage plots of the duplicates to see if they are evenly distributed or represent only a few genes. After aligning to the genome with HISAT2, I've flagged the duplicates with Picard and used Samtools depth to get the counts (of non-duplicates). I then plot them in R.
My problem is that I cannot get counts for just the duplicate reads. Here is what I've tried so far:
Extract duplicates from Picard MarkDuplicates output.
samtools view -h -f 0x400 input_markdups.bam > dups.sam
Calculate counts. Produces empty file with no errors.
samtools depth dups.sam > dups_depth.txt
Thinking that depth may be ignoring the duplicates (i.e. all reads), I removed the flags.
picard RevertSam I=dups.sam O=dups_revsam.sam REMOVE_DUPLICATE_INFORMATION=true
Once again, samtools depth runs for a minute or two and then ends with a blank file and no errors.
What am I missing? Is there an easier way to get counts for duplicates? Any help will be greatly appreciated.
Note: This is my first post here, so please let me know if I have deviated from standard procedures. Thank you!
