Hello
I have BAM-full file with reads mapped to "human and mouse" chromosome file. Now I would like to extract reads mapped only to "mouse" (means not mapped to human chromosome". This is the protocol I am using:
- From BAM-full, extract reads mapped to Mouse chromosome and convert to fastq1
- From BAM-full, extract reads mapped to Human chromosome and convert to fastq2
- Extract read ids from both fastq1 and fastq2 and identified non-overlapped "reads" ids from fastq1
I would like to know if there is a samtools command or software to apply these filtration steps?
Thanks a lot in advance
2 answers
I would recommend doing the alignment directly with bbsplit from BBTools with ambiguous2=toss and you are done.
While using bbsplit is a painless way of doing this once, if you do want to follow your original line of thought you can do the following.
Split human/mouse reads into separate BAMs using ideas from:
How To Split A Bam File By Chromosome
Extract ONLY chromosomes 1-22 from bam file - removing extraneous chr annotations Hopefully you have unique chromosome names from mouse and human otherwise nothing is going to work.Once you have the BAMs. you can
samtools collate | samtools fastqto get fastq-format reads. See the options for those commands to deal with secondary mappings etc. Step 1 and 2 can possibly be combined via a single pipeExtract fastq headers by something line
grep "@M01923:976:00000000" file.fq(usezgrepif your files are compressed). Remove@at beginning of the header at the same time. example belowzgrep @first_part_of_headers file_R1_001.fastq.gz | sed 's/^@//g' > fastq_headers_R1Use
filterbyname.shfrom BBMap suite to find reads that are present in one or other file.
Log in to answer this question.