Hello,
I am attempting to align sequence reads (fastq) to the PhiX genome in order to remove reads that map to the PhiX genome before going into my assembly. I have performed the alignment using BWA, and then extracted the unmapped reads from the alignment by using this command:
samtools view -f4 aln_B51_2phix.sam > aln_B51_2phix.unmapped.sam
I then convert the reads to fastq using this command:
cat aln_B51_2phix.unmapped.sam | grep -v ^@ | awk 'NR%2==1 {print "@"$1"_1\n"$10"\n+\n"$11}' > unmapped/aln_B51_2phix_0.fastq
and
cat aln_B51_2phix.unmapped.sam | grep -v ^@ | awk 'NR%2==0 {print "@"$1"_1\n"$10"\n+\n"$11}' > unmapped/aln_B51_2phix_0.fastq
I have a script which typically adds /1 /2 to the headers, and I added these to the reads before I did the alignment, but after the alignment they are no longer included in the fastq headers. I am unable to add the /1 /2 headers to these fastq files, even though they should be in the same format as the input fastq files for the alignment at this point.
I need the headers to include the /1 /2 designation for forward and reverse reads in order to perform the assembly. Does anyone know if I am doing something incorrectly or if there is a much simpler way to do all this?
Thanks!
•
written
by
p.migeon •
2 answers
Perhaps this might help you? But check the code first, since I did not test it.
samtools view -S -f4 aln_B51_2phix.sam | perl -F'\t' -ane 'chomp($F[10]); $mate=($F[1]&64)?"/1":"/2"; print "\@$F[0]$mate\n$F[9]\n+\n$F[10]\n";' | gzip -c > unmapped/aln_B51_2phix_0.fastq.gz
Ok, I didn't know that. How can an unmapped read be on the reverse strand? And what mappers do give me these reversed unmapped reads in the bam?
Hm... you're right. If one mate is mapped and the other mate not. Never thought about that issue. Thanks for clarifying it.
It should be noted that this is (finally) changing, see here.
Might also be useful to know that samtools has an option to generate fastq reads from BAM files.
Usage: samtools bam2fq [-nO] [-s <outSE.fq>] <in.bam>
So you can always:
samtools view -f 4 -Sb INPUT.sam > UNMAPPED.bam
samtools bam2fq -s Single_end.fastq UNMAPPED.bam > Paired_end.fastq
This should generate your Paired_end.fastq file with paired end reads with /1 and /2 appended. And if it so happens that you have single end reads, it will output them in Single_end.fastq
•
written
by
Adrian Pelin
cool, I didn't know that sub-command.
I could swear it wasn't there a while ago, but I am glad they implemented it, I really like how it works.
Nice one!
But it should be doable in one call:
samtools view -f4 -b INPUT.bam | samtools bam2fq - > Paired_end.fastq
Makes sense, the magic of command line! Thanks for the tip.
Log in to answer this question.
-
snakemake error: MissingInputException
written by tomilova.ekaterina.rp •
-
Snakemake. Problem with awk
written by kamanovae •
-
How to extract reads with exact match with BBDuk?
written by camelest •
-
Multiple commands in a for loop
written by lamia_203
-
Command for counting gene model after each round of MAKER
written by Elizabeth •
-
awk command in shell script
written by mittu1602
-
How to confirm forward and reverse primer matches the beggining of each read and keep only reads wi…
written by dzisis1986 •
-
Forum: Fastest way from SRA to fastq PE reads
written by Avi Srivastava
-
Extract read count for specific allele from a GATK VCF file
written by Sandeep
-
linux command in a perl script problem (sam to fastq)
written by rob234king
As a side note, BBMap can output fastq streams directly, making the process easier:
Though I generally use BBDuk for phiX removal instead, as it's faster:
In each case the original read names are maintained. You can also simplify arguments like this:
...as long as the files are named 1 and 2.