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$
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$ – user1738234 Jun 8 '20 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 '20 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 '20 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$ – user1738234 Jun 8 '20 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$ – thomas duge de bernonville Jun 8 '20 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$ – thomas duge de bernonville Jun 9 '20 at 5:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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