Unmapped reads to fastq (hisat2 / samtools / bedtools) not matching pair-end reads
I want to obtain a fastq file containing reads (RNASeq) that failed to map to a reference genome.
hisat2 -p 3 -x index -1 TD1.R1.fastq.gz -2 TD1.R2.fastq.gz -S out.sam --un-conc out.un.sam --summary-file out.info
I obtained that way 4 files (.sam + .info + out.un.1.sam + out.un.2.sam). Those ".un" files seems like the paired end reads that failed to map. The problem is that the bam2fastq scripts out there want only one bam and create two fastq files.
$ bedtools bamtofastq -i aln.qsort.bam \
-fq aln.end1.fq \
-fq2 aln.end2.fq
What can I use to take both .un files and create the two fastq files that I need? thanks!
Update: At the time I do the following
#align against first genome
hisat2 -p 3 -x index/svevo -1 reads/TD1.R1.fastq.gz -2 reads/TD1.R2.fastq.gz -S data/landic_svevo_1.sam --un-conc data/landic_svevo_1.unmapped.sam --summary-file data/landic_svevo_1.sam.info
# unmapped files by hisat2 are broken, error when converting to bam
#get unmapped reads with samtools
samtools view -h -f 4 data/landic_svevo_1.sam > data/landic_svevo_1.un.sam
#convert to bam
samtools view -b -S data/landic_svevo_1.un.sam > data/landic_svevo_1.un.bam
#get fastq from unmapped bam
samtools fastq -1 data/landic_svevo_1.un.1.fastq -2 data/landic_svevo_1.un.2.fastq data/landic_svevo_1.un.bam
#align unmapped against second genome
hisat2 -p 3 -x index/zavitan -1 data/landic_svevo_1.un.1.fastq -2 data/landic_svevo_1.un.2.fastq -S data/landic_un_zavitan_1.sam --un data/landic_un_zavitan_1.un.sam --summary-file data/landic_un_zavitan_1.sam.info
But get the next error with hisat2 at the end:
Error, fewer reads in file specified with -1 than in file specified
with -2 terminate called after throwing an instance of 'int' Aborted
(core dumped)
• 634 views
•
link
updated
by
Ram
•
written
by
juan.crescente •
•
written
by
juan.crescente •
0 answers
No answers yet.
Log in to answer this question.
More posts like this
-
Read Classification from ONT
written by maevalefeuvre •
-
How to obtain unmapped reads from the paired end samples using hisat2
written by otieno43 •
-
Converting Bam file to Fasta (Zipped)
written by Eliveri •
-
unmapped reads from 10X visium spaceranger bam file
written by Omics data mining
-
extract unmapped reads from the paired end samples using hisat2
written by winmoon6 •
-
How to plot alignment statistics; how many reads are mapped to the genome?
written by anamaria
-
bowtie2 --un vs --un-conc
written by blur
-
how can I align several fastq files with hisat2 and get htseq-count
written by Learner
-
Problem with Saminfo files
written by chowdhury •
-
Bowtie2 --un-conc output empy
written by t2
You could use(orsamtools fastqreformat.shfrom BBMap suite) to convert the two.unfiles to fastq and then userepair.shfrom BBMap suite to sync them up, so you can separate any singleton reads and properly pair the rest.Edit:
samtools fastqneeds BAM files. I have a complete answer usingreformat.shbelow.Any ideas why the unmapped file gives this error?
SAM format file is already text. You don't need to use
samtools view. Useless/cat/moreif you want to look inside.You should not need to convert it to BAM to usesamtools fastq, if that is what you intended to do.seems like samtools fastq only wants bam files, I've updated the question
Have you tried
samtools bam2fqfor extracting reads from.bamfile?http://www.htslib.org/doc/samtools-1.1.html
You are right. Use
reformat.shfrom BBMap suite in that case.That is because your
read 1/read 2files are no longer in sync.Use
repair.shfrom BBMap suite to fix that issue.seems like the output of hisat2 --un-conc is broken.
Try this out to see if it works: C: Unmapped reads to fastq (hisat2 / samtools / bedtools) not matching pair-end rea