{"id":57,"date":"2015-08-03T16:24:40","date_gmt":"2015-08-03T14:24:40","guid":{"rendered":"https:\/\/blogs.fu-berlin.de\/gruenstaeudl\/?p=57"},"modified":"2015-08-25T17:36:27","modified_gmt":"2015-08-25T15:36:27","slug":"filtering-out-unpaired-raw-reads-from-illumina-data","status":"publish","type":"post","link":"https:\/\/blogs.fu-berlin.de\/gruenstaeudl\/2015\/08\/03\/filtering-out-unpaired-raw-reads-from-illumina-data\/","title":{"rendered":"Filtering out unpaired raw reads from Illumina data"},"content":{"rendered":"

Working at the intersection<\/strong><\/p>\n

I have always wondered why an Illumina machine occasionally generates “unpaired” paired-end reads (i.e., when you receive an R1 but no corresponding R2, or vice versa). While I am waiting for a satisfactory answer, I would like to remove any unpaired reads from my data set in the meantime.<\/p>\n

Superficially, this aim may appear to be achieved through joining paired-end Illumina reads, which can be done through several software tools (e.g., FastqJoin<\/a>). In reality, however, all I wish to do is to filter out unpaired paired-end reads. Such a task is best completed by some shell scripting.<\/p>\n

First, I define my infiles.<\/p>\n

INF1=infile_R1.fastq
\nINF2=infile_R2.fastq
\n<\/code><\/p>\n

Second, I check if any unpaired paired reads are present.<\/p>\n

c1=$(wc -l $INF1 | awk '{print $1}')
\nc2=$(wc -l $INF2 | awk '{print $1}')
\nif (($c1 == $c2)); then echo \"TRUE\"; else echo \"FALSE\"; fi
\n<\/code><\/p>\n

Third, I extract the sequence names (but not the actual sequences, nor the quality scores).<\/p>\n

cat $INF1 | grep '^@' | grep -v ',' > $INF1.headers
\ncat $INF2 | grep '^@' | grep -v ',' > $INF2.headers
\n<\/code><\/p>\n

Fourth, I remove the unique pair-identifiers from the end of all sequence names. (This may vary depending on your data.)<\/p>\n

sed -i 's\/\\_1\\:N\\:0\\:3\/\/g' *.headers
\nsed -i 's\/\\_2\\:N\\:0\\:3\/\/g' *.headers
\n<\/code><\/p>\n

Fifth, I list all intersecting elements.<\/p>\n

comm -12 $INF1.headers $INF2.headers > intersect.headers
\n<\/code><\/p>\n

Sixth, I extract those reads from the input files, whose sequence names appear in the intersection. (Unfortunately, this step is slow. If anyone can come up with faster code, please write me.)<\/p>\n

for i in $(cat intersect.headers); do cat $INF1 | grep -A3 ^$i >> $INF1.intersect; done &
\nfor i in $(cat intersect.headers); do cat $INF2 | grep -A3 ^$i >> $INF2.intersect; done &<\/code><\/p>\n

Seventh, I perform some file hygiene.<\/p>\n

rm *.headers<\/code><\/p>\n

You then continue to work with the filtered files $INF1.intersect<\/em> and $INF2.intersect<\/em>, respectively.<\/p>\n

 <\/p>\n

EDIT:<\/em> Upon bouncing some ideas off of Stackoverflow<\/a>, I found a much faster code to conduct step 6 of the above pipeline (i.e., extracting those reads from the input files, whose sequence names appear in the intersection). Hence, the following code should be used instead:<\/p>\n


\ngrep -A3 -Ff intersect.headers $INF1 >> $INF1.intersect &
\ngrep -A3 -Ff intersect.headers $INF2 >> $INF2.intersect &
\n<\/code><\/p>\n

 <\/p>\n

 <\/p>\n","protected":false},"excerpt":{"rendered":"

Working at the intersection I have always wondered why an Illumina machine occasionally generates “unpaired” paired-end reads (i.e., when you receive an R1 but no corresponding R2, or vice versa). While I am waiting for a satisfactory answer, I would like to remove any unpaired reads from my data set in the meantime. Superficially, this […]<\/p>\n","protected":false},"author":2306,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[57598],"tags":[],"class_list":["post-57","post","type-post","status-publish","format-standard","hentry","category-bioinformatics"],"_links":{"self":[{"href":"https:\/\/blogs.fu-berlin.de\/gruenstaeudl\/wp-json\/wp\/v2\/posts\/57","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/blogs.fu-berlin.de\/gruenstaeudl\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.fu-berlin.de\/gruenstaeudl\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.fu-berlin.de\/gruenstaeudl\/wp-json\/wp\/v2\/users\/2306"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.fu-berlin.de\/gruenstaeudl\/wp-json\/wp\/v2\/comments?post=57"}],"version-history":[{"count":10,"href":"https:\/\/blogs.fu-berlin.de\/gruenstaeudl\/wp-json\/wp\/v2\/posts\/57\/revisions"}],"predecessor-version":[{"id":72,"href":"https:\/\/blogs.fu-berlin.de\/gruenstaeudl\/wp-json\/wp\/v2\/posts\/57\/revisions\/72"}],"wp:attachment":[{"href":"https:\/\/blogs.fu-berlin.de\/gruenstaeudl\/wp-json\/wp\/v2\/media?parent=57"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.fu-berlin.de\/gruenstaeudl\/wp-json\/wp\/v2\/categories?post=57"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.fu-berlin.de\/gruenstaeudl\/wp-json\/wp\/v2\/tags?post=57"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}