2
$\begingroup$

I'd like to fuse ~50 pairs of corresponding fastq files on per sequence bases. I'd be happy for a solution as a Linux script, R script (if feasible with huge files), or any free software that runs preferably under Debian Linux (otherwise Windows). Aim is preprocessing for a Linux pipeline for NGS data. Ideally both the sequence and quality data are fused.

The aim of this request is the following: there is an existing pipeline for downstream mapping of immunoreceptor sequences (migec/mixcr/vdjtools etc.). Previously we got Illumina reads (read 1, read 2 = reverse), with the index and UMI sequences attached to one of the reads => can be alligned and processed to habe 1 biological useful sequence with UMI and index at one end. NOW we get read 1, read 2 and seperately "read 3" which is just the UMI+Index reads, and which was attached by the core facility to read 1 using software we do not have access to. Since the previous pipeline works well, we prefere to use it in the future, but are open for suggestions how to use the index reads with migec/mixcr/vdjtools otherwise.

Each of the 2 files contains the same number (millions) of short (~300 bp) and corresponding DNA reads:

File_A.fastq:

  • @record1 / aaa... / +... / ...
  • @record2 / ...
  • @record3 / ...
  • ... (millions)

File_B.fastq:

  • @record1 / ttt... / +... / ...
  • @record2 / ...
  • @record3 / ...
  • ... (millions)

The target file should look like:

  • @record1 / aaa...ttt... / +... / ...
  • @record2 / ...
  • @record3 / ...
  • ... (millions)

(of course the aaa...ttt... is just an example)

I was looking for a solution on Stackexchange and Google, but the solutions, I found do not work for me because

  • manual processing is not an option given the number of sequences
  • I need to fuse the files on a per sequence basis, i.e. a simple Linux command like cat file1 file2 is not what I need (as far as I understand, this would give me (record1_A record2_A.... record1B record2_B).
  • R, package shortread, readFastq + writeFastq would be perfect, but cannot handle such large numbers of input lines

Maybe I did not use the right search terms, but in any case I'd be happy for a solution

$\endgroup$
4
  • $\begingroup$ If you want to join the forward and the reverse sequence, why do you want one sequence to be after the other one? Simply it isn't that way in the cell so most programs don't aim to do that. If these are not forward and reverse sequence, did you try perl, python or awk? Are you open to a solution in those languages? $\endgroup$
    – llrs
    Aug 15, 2018 at 15:56
  • $\begingroup$ Thanks - I added an explanation as second paragraph. These are of course not fw/rev reads. I have unfortunately not used python, perl or awk so far (but am open if it comes with an easy explanation) $\endgroup$
    – Martin
    Aug 15, 2018 at 17:02
  • $\begingroup$ From my understanding of the problem (which might well be wrong): If the sequences are in the same order (i.e. lines correspond to each other) one very simple approach would be just pasting the files together with paste -d "" File_A.fastq File_B.fastq > target.fastq. This would of course require some clean up first. $\endgroup$ Aug 16, 2018 at 9:01
  • $\begingroup$ @SebastianMüller: Yes, seq are corresponding and in the same order. And they have the same correspnding names in both files (sorry - I made this more precise in the question), i.e. only the sequence and quality data (i.e .lines 2 + 2n), but not the names (lines 1 + 2n) should be fused. Is there a way to do that, or afterward replace fused names with the original ones by another command? (I am new to Linux, was looking for sth. like this but could not find). If sth. like this works with 10+ mio seq per file, this would be the desired answer! $\endgroup$
    – Martin
    Aug 20, 2018 at 22:15

2 Answers 2

1
$\begingroup$

Based on the comments of the question I thought to put together this answer using GNU coreutils which should exist on pretty much all systems.

The idea is to paste the 2 files together with paste (assuming they are in the same order and each line correspond to the same line in the other file).

Since only the sequence and quality of both files should be kept, the lines starting with @ and + need to be made empty for one of the files. This can be achived with sed, e.g. sed 's/^@.*//' empties all lines starting with @ but keeps the empty line.

To modify the second file on the fly, one can using process substitution <(). Everything strung together should be something like this:

file1=File_A.fastq
file2=File_B.fastq
paste -d "" $file1 <(cat $file2 | sed 's/^@.*//' | sed 's/^+.*//') > target.fastq

I haven't tested it thoroughly, but could be a good starting point.

If you have a lot of files to paste, you might however want to process them first (not on the fly) using:

for file in $(ls *.fastq); do 
  cat $file | sed 's/^@.*//' | sed 's/^+.*//' > ${file%.fastq}_cleaned.fastq
done

This creates File_A_mod.fastq from File_A.fastq etc. Best to move them in a separate directory (make sure the first one is the original instead of the _cleaned one since it serves as the base). You should have the following in the folder: File_A.fastq File_B_cleaned.fastq File_C_cleaned.fastq .... Those can be than pasted with paste -d "" *.fastq > target.fastq.

$\endgroup$
0
$\begingroup$

What you actually want to do is overlap the reads, for which you can use tools like FLASH or BBMerge (I've had better luck with FLASH).

$\endgroup$
6
  • $\begingroup$ Thanks, the 2 sequences I want to process are NOT fw/rev or any other overlapping sequences (see new explanation above) $\endgroup$
    – Martin
    Aug 15, 2018 at 17:04
  • $\begingroup$ Ah, but at least mixcr will accept the two fastq files as they already are. $\endgroup$
    – Devon Ryan
    Aug 15, 2018 at 17:10
  • $\begingroup$ Yes, but somehow I need to keep the UMI/index from the 3rd file... If there is a way to do this with mixcr&co, that would be perfect, but I am not aware of that possibility. $\endgroup$
    – Martin
    Aug 15, 2018 at 17:41
  • $\begingroup$ The samples should have been split according to index (I assume this is an Illumina sample index) and the UMI is typically added to the read name (I or someone else can produce a few lines of code to do that for you). $\endgroup$
    – Devon Ryan
    Aug 15, 2018 at 17:42
  • $\begingroup$ Thanks a lot. Yes, it is Ilumina. Unfortunately sequences are not split for index, and index +UMI are not included in filenames (will add examples as soon as I get faster internet), but Migec can list each index in a barcode file and demultiplex according to index, if supplied with 2 (not 3) files, with UMI+index added to e.g. read 1... (used to work perfectly with this pipeliine as long as we had index+UMI added to read 1)... $\endgroup$
    – Martin
    Aug 16, 2018 at 8:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.