Hello communauty, the question is in the title: How to Trim bam file reads to remove all nucleotide over a certain length.
For fastq it sounds possible with cutadapt, but I need to do that on a bam file.
On this post Trim bam file reads to the same length they talk about reformat.sh but it just remove read over/under a certain length.
I have tried samtools ampliconclip thinking I could cheat to perform what I want wihtout success...
I have tried Cycle based clipping from GATK but it remove all reads (WellformedReadFilter remove everything...)
Any idea?
2 answers
Just as a thought: That (trimming sequences) will completely mess up the synchronization between the CIGAR string and the sequences. I could well imagine that downstream tools might refuse to use these "corrupted" files, maybe even samtools as I think it runs some basic validation routines that might choke here, not sure about that, you'll see :)
And of course reads that run reverse will have their fronts chopped off, and this will mean their mapping position is incorrect, and if this is paired-end data, that means the mate data of the mate will have to be fixed as well. And the mapping quality scores will be wrong, since they were calculated with longer reads.
I hope there is no clipping going on.
It's far easier to document what you are doing if you use established tools to do established things, then to make your own tools and have to work out all the bugs and all the corner cases. There are tested tools to regenerate the fastqs, and tested tools to trim.
I agree it is whay I was looking for a known tool alredy performing this task. I'm surprised that it does not exist.
To be clear what I want to perform: All my reads start at the same position, same strand. I want an even coverage all along the investigated region so I want to trim the ends to shorten the region to a specific position.
----------------
------------
-------------
should become
------------
------------
------------
Yes it does.
$ samtools view new.bam | head -2
NS500177:19:H2HLYXXXX:1:11104:17976:3958 97 chr1 39850 13 5=1X4=1X1=1X18=2X1=1X27=1X3=91D1X7=2X chr22 26122134 0 GAGCCGAGATTGCGCCACTGCACTCCAGCCTGGGCGACAGAGCGAGACTCTGTCTCAAAAAAAAAAAGGAGGGGGG <<AAAFFF<FFFFAFF7FFAFAFFFFFFFFFFFFFFFFFF<FF..FAFFFAFFFFFFFFFFFFFFF.F<FFFFFFF NM:i:10 AM:i:13 XS:A:+
NS500177:19:H2HLYXXXX:1:11103:20522:14693 129 chr1 61792 23 6=72D1=1X22=1X8=1I8=1X27= chr22 26121634 0 CCCTCCTTTTTTTTTTTGAGACAGAGTCTCGCTCTGTCGCCCAGGCTGGAGTGCAGTGGCGCAATCTCGGCTCACT AAAAAFFFFFFFFFF<FF<F.FFFFFFF<FA<<FFFAF7FFFFF.FFFFFFFF7AFFFFFFFFF..)FFFAFF.FF NM:i:4 AM:i:23 XS:A:-
$ reformat.sh -Xmx5g in=new.bam out=trim.bam forcetrimright=50 ref=genome.fa
$ samtools view trim.bam | head -2
NS500177:19:H2HLYXXXX:1:11104:17976:3958 97 chr1 39850 13 5=1X4=1X1=1X18=2X1=1X16= chr22 26122134 0 GAGCCGAGATTGCGCCACTGCACTCCAGCCTGGGCGACAGAGCGAGACTCT <<AAAFFF<FFFFAFF7FFAFAFFFFFFFFFFFFFFFFFF<FF..FAFFFA XS:A:+
NS500177:19:H2HLYXXXX:1:11103:20522:14693 129 chr1 61792 23 6=72D1=1X22=1X8=1I8=1X2= chr22 26121634 0 CCCTCCTTTTTTTTTTTGAGACAGAGTCTCGCTCTGTCGCCCAGGCTGGAG AAAAAFFFFFFFFFF<FF<F.FFFFFFF<FA<<FFFAF7FFFFF.FFFFFF XS:A:-
Log in to answer this question.
-
Samtools ampliconclip hard clip not giving same results as soft clip
written by Kevin •
-
Extract gRNA sequence using cutadapt
written by Swimming bird •
-
How to calculate z-Score over bins in bam file
written by conradinbaumi •
-
Removing adaptors from smallRNA seq data
written by juan.crescente •
-
Incorrect cycle number entered on sample sheet
written by harks •
-
rnaseq sequence analysis and trimming
written by Björn •
-
biopython script to preprocess raw RNA-seq reads (quality filtering, polyA and adapter trimming)
written by shanasabri •
-
Many 0 length reads after trimming
written by marcon •
-
About The Priority Of Trimming
written by Tonyzeng •
-
Bam Files And Reads Quality Check
written by madkitty
@Pierre will likely have something either pre-written or ... Did you check his
jvarkit(LINK) page?