2
$\begingroup$

I have a FASTA file that has reads like this:

>SRR5655563.745 745 length=126 (type=T,start=2,end=12,length=11,identity=81.8182%) CttttatttgttGTGTGTGAGGTTTGATTTGATGGGAAAATATCTTGAATCTGCGGCGAGGTTGGAAGAGCTATCGCGGATTGTGTCATCTGCTGCGAAGCCCAATAGGTCAAAGGGAACGCTACC

As you can see, there is information about the start "start=2" and end "end=12" of a sequence within the header. I would like to slice the sequence like [start:end] and keep the rest of it. e.g

The part I would like to trim:

Ctttggtttcctttt

And after trimming, I would like to keep the rest of the read:

GTGTGTGAGGTTTGATTTGATGGGAAAATATCTTGAATCTGCGGCGAGGTTGGAAGAGCTATCGCGGATTGTGTCATCTGCTGCGAAGCCCAATAGGTCAAAGGGAACGCTACC

So the new FASTA file will have reads like this:

>SRR5655563.745 745 length=126 (type=T,start=2,end=12,length=11,identity=81.8182%) GTGTGTGAGGTTTGATTTGATGGGAAAATATCTTGAATCTGCGGCGAGGTTGGAAGAGCTATCGCGGATTGTGTCATCTGCTGCGAAGCCCAATAGGTCAAAGGGAACGCTACC

I have written a python script that takes a FASTA file as an input and outputs a trimmed FASTA file, like this:

tailtrim.py --fasta input.fasta --output output.tail.trimmed.fasta

The code does trim the reads where I want it to.

import numpy as np
import argparse
import re

def main():
    parser = argparse.ArgumentParser()

    parser.add_argument("--fasta", type=str, required=True, help="input fasta file")
    parser.add_argument("--output", type=str, required=True, help="output tail.trimmed.fasta file")

    args = parser.parse_args()

    out_fasta = open(args.tail.trimmed.fasta, 'w')

    for line in open(args.fasta):
        header = line.rstrip()
        read = next(args.fasta).rstrip()

    #Extract header into a list
        header_lst = []
        header_lst= header.split(" ")
    
        tmp_start_lst = []
        tmp_end_lst = []
        tmp_lst = []

        for i in header_lst:
            if i.startswith("(") and i.endswith(")"):
                tmp_lst.append(i.split(","))
        #print(tmp_lst)

    # Extract the start and end numbers in a list
        for i in tmp_lst:
            for str in i: 
                if str.startswith("start="):
                    tmp_start_lst.append(str)
                if str.startswith("end="):
                    tmp_end_lst.append(str)

        pattern = re.compile(r'\d+$')

        start_lst = []
        for str in tmp_start_lst:
            match = re.search(pattern, str)
            start_lst.append(int(match.group()))

        end_lst = []
        for str in tmp_end_lst:
            match = re.search(pattern, str)
            end_lst.append(int(match.group()))

        #Extract reads into a list to manipulate i.e. trim
        read_lst = []
        read_lst.append(read)

        new_lst = []
        
        for str in read_lst:
            for num_start in start_lst:
                for num_end in end_lst:
                    s = str.replace(str[(num_start-1):num_end],"")
                    new_lst.append(s)

I do get the trimmed reads, but what I'm having trouble figuring out is how to write these trimmed reads with their corresponding headers to a new FASTA file.

Also if there is any way that I could improve my code or employ better strategy to approach the problem, I would love to know.

Best wishes.

Edit: I tried running the solution by the super helpful user @Steve.

The command I ran:

python3 tailTrimming.py --fasta tailtrim.fasta --output testtrim.trim.fasta

The command line error:

Traceback (most recent call last):
  File "/home/x/Desktop/fastq/tailTrimming.py", line 61, in <module>
    main()
  File "/home/x/Desktop/fastq/tailTrimming.py", line 58, in main
    SeqIO.write(trim_reads(fasta), out, "fasta")
  File "/home/x/.local/lib/python3.10/site-packages/Bio/SeqIO/__init__.py", line 478, in write
    count = writer_class(fp).write_file(sequences)
  File "/home/x/.local/lib/python3.10/site-packages/Bio/SeqIO/Interfaces.py", line 209, in write_file
    count = self.write_records(records)
  File "/home/x/.local/lib/python3.10/site-packages/Bio/SeqIO/Interfaces.py", line 193, in write_records
    for record in records:
  File "/home/x/Desktop/fastq/tailTrimming.py", line 20, in trim_reads
    result = match.groupdict()
AttributeError: 'NoneType' object has no attribute 'groupdict'

This is the whole input file I supplied:

>SRR5655563.9676 9676 length=126 (type=T,start=1,end=11,length=11,identity=81.8182%)
tttctttctttAGTTCCAATAGTCTTGACCACCTTTTATCAACAGAAGCTTGGCCTTGTCCAGTTATATATTGTATTGTACCCTCAGATTTCTCGCACAAT
>SRR5655563.10087 10087 length=126 (type=T,start=2,end=16,length=15,identity=80%)
GtttgtttcattttttCCAAGGAAGAAATTCGACCTTCCAGCGATTTGTGGAGGTCATTCATCTGTTGAGTAAAGGATTGGTGCGCGACCTTGTAATCCTTGTGGAGGGTGTAG
>SRR5655563.13902 13902 length=126 (type=T,start=2,end=16,length=15,identity=80%)
CtttttcttcttgtttGATGTAAAAAATCTTCGGTATCTAAAACGTATGAAACCTTAGATGAACTGCAAGCTATCTTACACTCTAGACCATTGATAATTT
>SRR5655563.18489 18489 length=126 (type=T,start=2,end=12,length=11,identity=81.8182%)
CtttttcttcttCTTGGAAGTTTTGGATTCTGCAGGGACATTTTCTTTCTTTTCTCCATCTCCATCAGTCTCAACACCTGTCTTATCCTGAGCACCTTGAGACTCATCTTGACCATCTTTACTTT
>SRR5655563.20782 20782 length=126 (type=T,start=2,end=15,length=14,identity=71.4286%)
CttcattttagttttCAACTCCTCTACCAACTGGGATGGTGCATTCAACTGTAGAACTGTTTCACGAATCTCCTTGAGAACACTGAAGTCTCCCTCAGTTAACTTCTTTTTCAAATTATTGA
>SRR5655563.23972 23972 length=126 (type=T,start=2,end=29,length=28,identity=75%)
CttttttgttgttctcttctttcttctttGCATGGCTTGATGATGACTGATCCGATTGATTATTTCTTGATAGAATATCAATGTCATCACCAATTGTTGAGT
>SRR5655563.31739 31739 length=126 (type=T,start=1,end=10,length=10,identity=80%)
ttttcttcttCTTGCATGAAAGTGGGTTTGGACCCTTAGCTCTTTTTCTCTTGAACTGAGGTTT
>SRR5655563.34244 34244 length=126 (type=T,start=2,end=14,length=13,identity=76.9231%)
CtttcctttattttCTTGGAGGTTAGCCAGGATATTGGGATCAAAGCATGTAGCAAAACACACCTTGCTATTTCTTGCTCGTACAGATACATGGGATAGGA
>SRR5655563.38334 38334 length=126 (type=T,start=2,end=12,length=11,identity=81.8182%)
CtttattgttttCTCATCACCCATATTATATTTCTTAAGGGCAGATTCCAACTCTTCCATTGTGATGTACCCGCTCTTATCATTGTCAAAATATTCAAAGGCTTTATATAGGCGATCCTCTCTTTC

And this is the whole output I have gotten:

>SRR5655563.9676 9676 length=126 (type=T,start=1,end=11,length=11,identity=81.8182%)
AGTTCCAATAGTCTTGACCACCTTTTATCAACAGAAGCTTGGCCTTGTCCAGTTATATAT
TGTATTGTACCCTCAGATTTCTCGCACAAT

Just one read, that's all. Whereas I would like to keep all the sequences because all of them have this start, end attribute. The sequence being on two lines in the output is not something I did. It is like this in the output file. So it's no longer FASTA file, right?

Thanks in advance for any help!

Improve this question
$\endgroup$
4
  • 3
    $\begingroup$ I think there's a typo in your example. The part to trim does not exist in the input sequence. Also, if you want to keep the remaining sequence (i.e. "the rest of it"), shouldn't the output include the sequence before the starting position? Or would you like to extract the sequence after the end position? $\endgroup$
    – Steve
    Dec 24, 2022 at 13:17
  • 2
    $\begingroup$ Just to add in terms of performance Python is slow (loops and stuff). re is slow even for Python. Although Perl is not popular it will outperform the solutions here on a single CPU. If there are performance issues it's worth posting back - ultimately it depends on the size of the file and sensitivity to runtime. Huge files can simply fall over (i.e. never complete). However, I agree with Steve's solution: what Python is very cool at is splicing a string and that is super important to the solution. This is the bit that is really missing in the above code IMO $\endgroup$
    – M__ ♦
    Dec 24, 2022 at 16:45
  • $\begingroup$ You made mistake describing the problem (sequence does not match 2-12 position intervals). If you update it, maybe you get more help. Also, it's not clear what kind of header you want to print in the record, exact same FASTA header? $\endgroup$
    – Supertech
    Dec 24, 2022 at 20:36
  • $\begingroup$ i have the same problem but looks like i found some interesting solution @Steve. However, i have want to ask... what does this mean? pattern = ( r'type=(?P<type>\w+),' r'start=(?P<start>\d+),' r'end=(?P<end>\d+),' r'length=(?P<length>\d+),' r'identity=(?P<identity>\d+(?:\.\d+)?)%' ) r'type= ??? what should i put here? r'start = ??? refers to the forward primer? r'end = reverse primer? r'length = the length of the final sequence after the trimming? what happens if there are sequences with different length afterwards? r'identity = the similarity between the primer and its sequence? $\endgroup$
    – Damianos
    Feb 10 at 21:19

1 Answer 1

Reset to default
3
$\begingroup$

Here's one way using Biopython and the SeqIO interface to read and write SeqRecord objects. Using a generator, we can produce (yield) trimmed sequence records that can then be written to the output handle. This avoids storing an entire list of trimmed sequences in memory and is usually the way to go, assuming you're dealing with large files containing thousands of records:

import argparse
import re

from Bio import SeqIO

def trim_reads(fileobj):

    pattern = (
        r'type=(?P<type>\w+),'
        r'start=(?P<start>\d+),'
        r'end=(?P<end>\d+),'
        r'length=(?P<length>\d+),'
        r'identity=(?P<identity>\d+(?:\.\d+)?)%'
    )

    for seq_record in SeqIO.parse(fileobj, "fasta"):
        match = re.search(pattern, seq_record.description)
        result = match.groupdict()

        start = int(result['start']) - 1
        end = int(result['end'])

        head = seq_record.seq[:start]
        tail = seq_record.seq[end:]

        seq_record.seq = head + tail

        yield seq_record

def get_argument_parser():

    parser = argparse.ArgumentParser()

    parser.add_argument(
        "--fasta",
         type=str,
         required=True,
         help="input fasta file",
    )
    parser.add_argument(
        "--output",
        type=str,
        required=True,
        help="output tail.trimmed.fasta file",
    )

    return parser

def main():

    parser = get_argument_parser()
    args = parser.parse_args()

    with open(args.fasta) as fasta, open(args.output, "w") as out:
        SeqIO.write(trim_reads(fasta), out, "fasta-2line")

if __name__ == "__main__":
    main()
Improve this answer
$\endgroup$
8
  • 1
    $\begingroup$ Its worth point out that re is notoriously slow and over large files parallelisation is needed across the Biopython object. $\endgroup$
    – M__ ♦
    Dec 24, 2022 at 15:21
  • 1
    $\begingroup$ No worries at all, @pubsurfted. Yes - re.search will return None if the search pattern could not be found. You'll need to add some code to handle these cases. You might need a different search pattern, or maybe you just want to ignore these records? If, in your question above, you can provide a few examples of where the search pattern fails, I might be able to suggest a solution. $\endgroup$
    – Steve
    Dec 27, 2022 at 11:18
  • 1
    $\begingroup$ @pubsurfted I've updated the regex above to handle the fact that the identity can also be a plain integer, i.e. the decimal places are optional. Also, FASTA files are often line wrapped at 60 characters, but this entirely optional. These are sometimes called multiline FASTA files. There's probably a way to change the line wrapping using Biopython. $\endgroup$
    – Steve
    Dec 28, 2022 at 14:25
  • 1
    $\begingroup$ @pubsurfted There is indeed, using the "fasta-2line" format. $\endgroup$
    – Steve
    Dec 28, 2022 at 14:32
  • 1
    $\begingroup$ It ran. I had to update the biopython package from 1.67 to 1.80. $\endgroup$
    – pubsurfted
    Dec 29, 2022 at 14:16

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.