2
$\begingroup$

I have a FASTQ file:

@NS500455:80:HG7TNBGXB:1:11101:17723:1055 1:N:0:ATCACG
ACTTANGTGTATGTAAACTTCCGACTTCAACTGTATAGGGATCCNAGCTCCAATTCGCCCTATAGTGAGTCGTAT
+
/AAAA#EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE#EEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
@NS500455:80:HG7TNBGXB:1:11101:8821:1057 1:N:0:ATCACG
ACTTANGTGATGTAAACTTCCGACTTCAACTGTATAATAAATATCTAGATCGGAAGAGCACACGTCCGAACTCCA

I want to remove the string ACTTAAGTGTATGTAAACTTCCGACTTCAACTG from it while retaining the sequences starting from TA onwards. I wrote grep "ACTTAAGTGTATGTAAACTTCCGACTTCAACTGTA" SRR_1.fastq | sed "s/ACTTAAGTGTATGTAAACTTCCGACTTCAACTGTA/TA/g", but it removes the headers of the FASTQ file as well.

The desired output would be

@NS500455:80:HG7TNBGXB:1:11101:17723:1055 1:N:0:ATCACG
TATAGGGATCCNAGCTCCAATTCGCCCTATAGTGAGTCGTAT
+
/AAAA#EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE#EEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
@NS500455:80:HG7TNBGXB:1:11101:8821:1057 1:N:0:ATCACG
TATAATAAATATCTAGATCGGAAGAGCACACGTCCGAACTCCA

How to do this. I dont want to use cutadapt.

Kindly help.

$\endgroup$

3 Answers 3

4
$\begingroup$

What application are you using that will be okay with a mismatching quality string?

Use trimmomatic or cutadapt, or bcl2fastq could probably do the job too. Those will also have the benefit of dealing properly with one-off errors, which awk and sed won't.

$\endgroup$
3
$\begingroup$

Skip the grep, use just the sed. Unless the header contains the long sequence as well (which I doubt it will), your command will work.

sed 's/^ACTTAAGTGTATGTAAACTTCCGACTTCAACTGTA/TA/' SRR_1.fastq > SRR_1.edited.fastq
$\endgroup$
3
  • $\begingroup$ It works but there are some sequences begining with: ACTTANGTGTATGTAAACTTCCGACTTCAACTG They are being retained. How to deal with them. Kindly guide. $\endgroup$
    – Amit
    Jun 8, 2020 at 15:18
  • $\begingroup$ That's because N does not match A,T,C or G as a literal. Replace all N with [ATCG] to remove these sequences. $\endgroup$
    – Ram RS
    Jun 8, 2020 at 17:10
  • $\begingroup$ It is worth considering thomas' point. If you'd like your QUAL lines to reflect the change to the SEQ lines. awk or bioawk might be better at this point. $\endgroup$
    – Ram RS
    Jun 9, 2020 at 14:04
2
$\begingroup$

If you need the quality vector for each read, then your code will create a shift because it does not cut the quality vector. You should rather use a trimming tool a systematically cut all 33 first characters. fastp will allow you to do it easily with the -f, --trim_front1 argument.

$\endgroup$
3
  • $\begingroup$ What should be the command line , I typed fastp -i SRR_1.fastq -o trimSRR1.fq -f 33 and it gives an error ERROR: trim_front1 (--trim_front1) should be 0 ~ 30, suggest 0 ~ 4. Kindly help. $\endgroup$
    – Amit
    Jun 8, 2020 at 15:50
  • $\begingroup$ strange! I did'nt know it was restricted to 30. You can run fastp twice with 20 in the first trimming and 13 in second one? not really convenient but it should work! $\endgroup$ Jun 8, 2020 at 15:59
  • $\begingroup$ you may alternatively had a look at Trimmomatic (usadellab.org/cms/?page=trimmomatic). It has a HEADCROP argument, but I don't know if it also restricted in length. $\endgroup$ Jun 9, 2020 at 5:14

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.