This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Union of unaligned fastq reads

Hi,

Raw data is a paired end fastq file.

Aligned it with Genome-1 (using STAR) and got unaligned R1 and Unaligned R2.

Also aligned it with Genome-2 (using STAR) and got unaligned R1 and Unaligned R2.

Please advise what would be the best method/tool to obtain the 'union' fastq of the unaligned reads from the unaligned of both the genome alignments.

Jeffin

rna-seq
• 2,733 views • link
written by Jeffin Rockey Jeffin Rockey

Do you really want the union or rather the intersection?

• 0 views • link
written by cschu181 cschu181

Union itself is the requirement

• 0 views • link
written by Jeffin Rockey Jeffin Rockey

Creating the union would be simply combining the unaligned files together? Just need to avoid duplicates.

• 0 views • link
written by WouterDeCoster WouterDeCoster

Assuming unaligned R1/R2 are .bamfiles:

samtools view unaligned.R1.genome1.bam | cut -f 1 > unaligned.R1.txt
samtools view unaligned.R1.genome2.bam | cut -f 1 >> unaligned.R1.txt
grep -A 3 -F -f <(sort -u unaligned.R1.txt) original.R1.fq | grep -v "\-\-" > union.unaligned.R1.fq

Analogously for .R2

• 0 views • link
written by cschu181 cschu181

Thanks.

Could you please provide some detail on whats happening in the grep line?

• 0 views • link
written by Jeffin Rockey Jeffin Rockey
  1. We fgrep (grep -F, i.e. take patterns from file) the original fastq file for all the read ids of the union of unaligned reads.
  2. -f takes the pattern file, which here (because I'm lazy) are given as a process substitution: <(sort -u unaligned.R1.txt) provides the output of the sort operation as a 'file' (the shell pretends this is a file). Alternatively, one could have done the sort -u into another file and used that filename as argument to the -f.
  3. grep -A 3 returns the line with the match + the following 3 lines. However, if there are multiple hits, then they will be separated by --, which is why the output of the first grep is piped into a grep -v, that ignores the -- lines.
• 0 views • link
written by cschu181 cschu181

Thanks a lot @cschu181 . Good to learn about the -A functionality in grep :)

• 0 views • link
written by Jeffin Rockey Jeffin Rockey

Hi, The best way to avoid duplicate fastq entries is the aspect what I am doubtful about.

If I write a small script it would be easily doable.But I wanted to know whether there is better method to combine while keeping duplicates away

• 0 views • link
written by Jeffin Rockey Jeffin Rockey

0 answers

No answers yet.

Log in to answer this question.

More posts like this