text:
Find open reading frames in a DNA sequence
I have a fast file containing DNA sequences. I've already cleaned the file and has only DNA sequences with no gaps or errors or white space characters. I want to find open reading frames in the fasta file that are longer than 30 bp.
Here's a little description: An open reading frame is defined by a sequence that

start with ATG,
ends with TAG, TAA, or TGA,
is a multiple of three.

For each open reading frame, I want to print the start position in the sequence, the length of the open reading frame, and the open reading frame sequence to the screen. Then I want to print the total number of open reading frames found. I would appreciate some help.
I can't use libraries.
Also, here's the DNA sequence I am working with:
GCATGCAATACAGTGACATATATATATACCCTAACACTACCCTAACCCTACCCTATTTCA
ACCCTTCCAACCTGTCTCTCAACTTACCCTCACATTACCCTACCTCTCCACTTGTTACCC
TGTCCCATTCAACCATACCACTCCCAACCACCATCCATCCCTCTACTTACTACCACCAAT
CAACCGTCCACCATAACCGTTACCCTCCAATTAGCCATATTCAACTTCACTACCACTTAC
CCTGCCATTACTCTACCATCCACCATCTGCTACTCACCATACTGTTGTTCTACCCTCCAT
ATTGAAACGTTAACAAATGATCGTAAATAATACACATATACTTACCCTACCACTTTATAC
CACCACACATCACATGCCATACTCACCTTCACTTGTATACTGATATGCCATACGCACACG
GATGCTACAGTATATACCACTCTCAAACTTACCCTACTCTCACATTCTACTCCACTCCAT
GACCCATCTCTCACTAAATCAGTACTAAATGCACCCACATCATTATGCACGGCACTTGCC
TCAGCGGTCTATACCCTGAGCCATTTACCCATAACTCCCACGATTATCCACATTTTAATA
TCTATATCTCATTCGGCGGGCCCAAATATTGTATAACTGCTCTTAATACATACGTTATAC
CACTTTTGCACCATATACTAACCACTCAAATTATATACACTTATGCCAATATAACCAAAA
AATTACCACTAAAATCACCTAAACATAAAAATTATTTATCTTTCAACTTTACGAAATAAA
CACACTCAATTGCGTATCTATACCACCATGACGTCATTAACGTAAAAGTTCCTTAATATT
ACCATTTGCTTGAACGGATACCATTTCAGAATATTTCTAACTTTCACAGACCATACATTA
GAATAATATGCCACCTCACTGTCGTAACACTCTATATTCACCGAGAAACAATACGGTAGT
GGCTCAAACTCATGCCGGTGCTATGATACAATTGTATCCTATTTCCATTCTCATATGCTA
TCCGCAATATCCTAAAAGCATAACTGATGCATCTTTAATCTTGTATGTGACACTACTCAT
ACGAAGAGACTATATCTAAAGAAGACGATACAGTGATATGTACGTTGTTTTTGTAGAATT
ATAATGAAACGTCAAATAACCCTACTATATTATAACTTATCAGCGCCGTATACTAAAACG
GACGTTACGATATTGTCTCACTTCATCTTACCACCCTCTATCTTATTGTTGATAAAACAC
TAACCCCTCAGCTTTATTTCTAGTTACAGTTACAACAAACTATCCCAAACCATAAATCTT
AATATTTTAGGTGTCAAAAAATGAGGATCTCCAAATGAGAGTTTGGTACCATGACTTGTA
ACTCCACTACCCTGATCTGCAATCTTGTTCTTAGAAGTGACGCATACTCTATATGGCCCG
ACGCGACGCGCCAAAAAATGAAAAAAGAAGCAGCGACTCATTTTTATGGAAGGACAAAGT
GCTGCGAAGTCATACGCTTCCAATTTCATTATTGTTTATTGGACATACTCTGTTAGCTTT
ATTACCGTCCACGCTTTTTCTACAATAGTGTAAAAGTTTCTTTCTTATGTTCATCGTATT
CATAAAATGCTTCACGAACACCGTCATTGATCAAATAGGTTTATAATATTAATATACATT
TATATAATCGGCGGTATTTATATCATCAAAAAAAGTAGTTTTTTATTTTATTTTTTCATT
ACTTTTCACTGTCTATGGATTTTCATTCGTAAAGGCATCACTCCCTAGTTTGCGATAGTG
TAGATACCGTCCTTGGATAGAGCACTGGAGATGGCTGGCTTTAATCTGCTGGAGTACCAT
GGAACACCGGTGATCATTCTGGTCACTTGGTCTGGGGCAATACCAGTCAACATGGTGGTG
AAGTCACCGTAGTTGAAAACGGCTTCAGCAACTTCGACTGGGTAGGTTTCAGTTGGGTGG
GCGGCTTGGAACATGTAGTATTGGGCCAAGTGAGCTCTGATATCAGAGACGTAGACACCC
AATTCCACCAAGTTGACTCTTTCGTCAGATTGAGCTAGAGTGGTGGTTGCAGAAAGCAGT
AGCAGCGATGGCAGCGACACCAGCAGCGATTGAAGTTAATTTGACCATTGTATTTGTTTT
GTTTGTTAGTGCTGATGTAATCTTAACAAGAAATAGTGAAATGAAAGCGCATACCTCAAA
GGCATATAGTTGAAGCAGCTCTATTTATACCCGTTCCTCCATCTGTCATCACTACTTAAA
CGATTCGTTAACAGACGCTCATTTAGCACCTCACATATTCTCCATATCTCATCTTTCACA
CAATCTCATTATCTCTATGGAGATGCTCTTGTTTCTGAACGAATCATACATCTTTCATAG
GTTTCGTATGTGGAGTATTGTTTTATGGCACTCATGTGTATTCGTATGCGCAGAATGTGG
GAATGCCTATTATAGGGGTGCCGGGGGGTGCCTTGAAAAACCCTTTTGCGCGCCTGTTAA
GTTTCCGTTTTCAGTCAAAAAGAATATCCGAATTTTAGATTTGGACCCTCGTTCAGAAGC
TTACTGTCTAAGCCATCATTTGGTGTGTCCTAAACGGTTTCCATGCAAAGCAACATCTTT
GTTACTCATTCCTGAAGGTTGAATAAAAAATAGCAACCATGTGCAGGGAGTCGTATATTG
TTAGATTCTAGAGACTTGTACGCATCATCAAAGCTGTAAATAGAATAAACATACGCAAGG
CGTCAAAAGTGCATAGTTAAGAAAATTCCTGACATGTGAAAATATGTGTTTATGAAATGT
GTCAAGGCCCGTCTATAGCGTAGTTAACCCCTCTGCAGGAGTAAGTGACTTTTTTTACGC
TCAAAAGGCAACGAGGGCACATACTTAAAAGTCATTTTCAAACACATCTGCAGTTTGCAA
CGACAGATAACAATATTATGATAGGATGGTATGATGTTATTAGCTGCCACATATTTTTAA
TAATAGCATTAGTCACGTCTCTTCAATTGTTGGGATGAAACTCTAAAATATCATTCCTTT
AGTAGTATTCCAGTTACCAGTATATTATCACATGCCGAAAAAGAAGATGACATAAAGATC
GACAAACAGTCTTCAAATATAATGGAAGCTGGAATGCAAGGATTGATAATGTAACAGGAT
ACTGAATGACAAAGTATAAATGAAAAAAAAAAAAAAAAAAAAGTAGTAATACTATTATGT
GGAAATACCGATTCCATTTTGAGGATTCCAATTGTTGGAATAAAAATCAACTATCATCTA
CTAACTAGTATTTACGTTACTAGTATATTATCATATACGGTGTTAGAAGATGACGCAAAT
GATGAGAAATAGTCATCTAAATTAGTGGAAGCTGAAACGCAAGGATTGATAATGTAATAG
GATCAATGAATATTAACATATAAAACGATGATAATAATATTTATAGAATTGTGTAGAATT
GCAGATTCCCTTTTATGGATTCCTAAATCCTGGAGGAGAACTTCTAGTATATCTACATAC
CTAATATTATAGCCTTAATCACAATGGAATCCCAACAATTACATCAAAATCCACATTCTC
TACACCAATACCATCGACGAGAGCTTCTAGTAAATTGTATACATAACAGTATAACCCTTA
CCAACAATGGAATCTCAAAGATTATTAAATTATTCACAGACTCTGAGGATTCGGGTAAAA
TAGGGTATTTAACTGGTTACCGGAAAGGTTTAGAAAATTCGTGGAGGGTTGGCCGAGTGG
TCTAAGGCGGCAGACTTAAGATCTGTTGGACGGTTGTCCGCGCGAGTTCGAACCTCGCAT
TCATTTTCCGGATCCACTTTATAGTGGGTCACTTCATCTCTTGTAGGGTTCACTGAAAGA
TAATCAATCAAATTATCATATGACCAATGGGTCACTTTTTTGATGAGACCCATGTTTATT
CTTCTATACGTTCGTGATACTGCACTTACGACTACCAGTAACATCAGAATAACCATAGAG
CTTATGCTATCTTGAAATATATCATCGCCGTAAAGTATGATGTGTAGCGCATTAGTCTTC
TGCATTCATATTATCAGTGGTATACCGTTTGGTATATGCTTTCCAATCATGAAGTGTTGG
CTCTACTCTATATAAATTCGACATCCATTGATTTTCTCAATAAATGAGTTACCCAAGTAA
GCTTTCCATTGATACAAGTGATCTACATTCTTGCGACGCCAAATTATAAAGCACTAAAAA
TCATATCATACCCAATTGCGGGCAGCACATGTATATATTATACACTACTTATAACTAACC
TTCTTCATATATAACAATATGCCTGATAATACTAATGCATTTAAATCATATCATGAAAGA
AAATTACATGGGTTTTATTGACATAATTGCATTTAGAATACATAAAATTCTAAAGAATTA
ATATATCCAAAGTATTAGACATAACCAAGAATAATAGTGAATAATTTTAGATTTTGTTAC
ATATAATTCTGCTTGCCTATCTCTTCCACTCTTTTCAAAACGTTGCATGTAAGCGTTACT
AATATTCCGCTTTATTTTGTTGCAATTCCTAATTTTTTCATTACATTATCTTGCGAGTAC
GGAAGCGATTAACGTTCTCCCAATAGAAGGAACAAACATAGATATTGAAGTTTTACTGCT
TTTGCTTACCTGACCTTTTTCAAATTTAATTTTTTCCCGCTAATAAGACCATAAACTACC
CCGAACCAAATTCTAAAAGATAGTCAGCTGGATTAGAGTTGTCATCTCCAAACATTAATT
TTGCATTATCTTCGGCTTCAATCAAATCGCCTGATAAGAACTCCTTTAATTCTTCATGAA
TGTTTGTATGTGGATGACTCTCCATAGTGCCAGCATGATTGTGGTTACCGACCGAATCAT
ATAACGGTGGCTCCCAATTGTGCAATTCAGCCTTACTATTTTGTTCATTCATCAAAGCAT
TTGGGACAGATCTAATATCTATAATTCTTTCCTCACTATTCTCGCTATTATTTTGCCCCG
AACTGGCATGGTGGTTATTGGTAAAAGGAGATAATGCTGCGACAGAACTTTTCTTCTCTT
CCAATTCTTTTATCTTCGCCAATAACGCCTTCTTTTTTCGTGCTTGTATTTCTAAAATTT
CGGCCAGGTATTGTAAATATTCGACCGTTCTATCCAAAATGATACCCTTATTTGGTTTGA
TTTGTTTACCTAGGTCATCGTAATTCAATAAAGATGGTGGAACCAACTGGCCGAGTTCTT
TTATCTTTTGCTTTATTAATTCTCTTCTTCTCCTTTCGACGGCATTATGAAACTCTCTTT
TGCGCCTCAGTTTCTCATCAGAAGTTAACCCGCCTAAAATCTTTGGAACACTTCCAGGTC
CTATATTTTCAGTCATATTGCTACTTATTGAAGTGTGTCTTGTTCTGGGTGTGTTTATGC
TACCATGCCTAAAAGAAGATGAAAGGAAACTTCCTGCACGAAATGACGATGATGGTGATC
TTACCTTTGGGGAATATGTGGAAGATACGGATGCTGGGCCCAAGCTTTGTGGATTATAAG
AAAATGATGATGAATATGTGTTTGGTGTCATCATATCAGAATTGATGCTAGAAGATAAAG
AGGAGCTCAAATCATCAGTTAAGTTATACATTGTGTCATCGGTGCCATGCTGAAATAAAA
AGTCATGTGGTAACTCAGCTTTTAAGCCTTCTTGTTGCTGTAATGGAGTATTCATTGCTC

"I can't use libraries." - why not? Is this an assignment question?

The answer looks like homework.

ok I have a Python attempt, used Biopython only to load your DNA sequence;
working only on your sequence ( not calculating the reverse complement) and as per your instruction an ORF is something in between an ATG and a stop codon, so I am using every ATG as a start even if preceded by another ATG, lenght of ORFs is >30bp meaning >90 bases. Code:
from Bio import SeqIO
# record = SeqIO.read("sample_test.fasta", "fasta")
# record = SeqIO.read("sample_test_2.fasta", "fasta")

record = SeqIO.read("sample.fasta", "fasta")

# print(record.seq)
# print(len(sequence))

sequence = record.seq

orfs = []

for i in range(len(sequence)-2):

    if sequence[i]+sequence[i+1]+sequence[i+2] == 'ATG':

        print('ATG at : ', i)

        orf = 'ATG'
        x = 3
        while sequence[i+x:i+3+x] not in ['TAG', 'TAA', 'TGA'] and i+x < len(sequence):
            orf += sequence[i+x:i+3+x]
            x += 3
        if len(orf) > 90:
            orfs.append((i, orf , len(orf)))

cnt = 1
for i in orfs:

    print('\nsequence number : ', cnt ,'\nstart : ', i[0] ,'\nframe : ',  ((len(sequence)-len(sequence[i[0]:]))%3)+1 ,
          '\nsequence : ', i[1],'\nlenght sequence :', i[2],'\n\n')
    cnt += 1

Result, start is with a sequence that goes from 0 to len(sequence)-1:
sequence number :  1
start :  2
frame :  3
sequence :  ATGCAATACAGTGACATATATATATACCCTAACACTACCCTAACCCTACCCTATTTCAACCCTTCCAACCTGTCTCTCAACTTACCCTCACATTACCCTACCTCTCCACTTGTTACCCTGTCCCATTCAACCATACCACTCCCAACCACCATCCATCCCTCTACTTACTACCACCAATCAACCGTCCACCATAACCGTTACCCTCCAATTAGCCATATTCAACTTCACTACCACTTACCCTGCCATTACTCTACCATCCACCATCTGCTACTCACCATACTGTTGTTCTACCCTCCATAT
lenght sequence : 300

sequence number :  2
start :  524
frame :  3
sequence :  ATGCACGGCACTTGCCTCAGCGGTCTATACCCTGAGCCATTTACCCATAACTCCCACGATTATCCACATTTTAATATCTATATCTCATTCGGCGGGCCCAAATATTGTATAACTGCTCTTAATACATACGTTATACCACTTTTGCACCATATACTAACCACTCAAATTATATACACTTATGCCAATATAACCAAAAAATTACCACTAAAATCACCTAAACATAAAAATTATTTATCTTTCAACTTTACGAAA
lenght sequence : 252

sequence number :  3
start :  907
frame :  2
sequence :  ATGCCACCTCACTGTCGTAACACTCTATATTCACCGAGAAACAATACGGTAGTGGCTCAAACTCATGCCGGTGCTATGATACAATTGTATCCTATTTCCATTCTCATATGCTATCCGCAATATCCTAAAAGCATAACTGATGCATCTTTAATCTTGTATGTGACACTACTCATACGAAGAGACTATATC
lenght sequence : 189

sequence number :  4
start :  982
frame :  2
sequence :  ATGATACAATTGTATCCTATTTCCATTCTCATATGCTATCCGCAATATCCTAAAAGCATAACTGATGCATCTTTAATCTTGTATGTGACACTACTCATACGAAGAGACTATATC
lenght sequence : 114

sequence number :  5
start :  1143
frame :  1
sequence :  ATGAAACGTCAAATAACCCTACTATATTATAACTTATCAGCGCCGTATACTAAAACGGACGTTACGATATTGTCTCACTTCATCTTACCACCCTCTATCTTATTGTTGATAAAACAC
lenght sequence : 117

sequence number :  6
start :  1370
frame :  3
sequence :  ATGACTTGTAACTCCACTACCCTGATCTGCAATCTTGTTCTTAGAAGTGACGCATACTCTATATGGCCCGACGCGACGCGCCAAAAAATGAAAAAAGAAGCAGCGACTCATTTTTATGGAAGGACAAAGTGCTGCGAAGTCATACGCTTCCAATTTCATTATTGTTTATTGGACATACTCTGT
lenght sequence : 183

sequence number :  7
start :  1432
frame :  2
sequence :  ATGGCCCGACGCGACGCGCCAAAAAATGAAAAAAGAAGCAGCGACTCATTTTTATGGAAGGACAAAGTGCTGCGAAGTCATACGCTTCCAATTTCATTATTGTTTATTGGACATACTCTGTTAGCTTTATTACCGTCCACGCTTTTTCTACAA
lenght sequence : 153

sequence number :  8
start :  1457
frame :  3
sequence :  ATGAAAAAAGAAGCAGCGACTCATTTTTATGGAAGGACAAAGTGCTGCGAAGTCATACGCTTCCAATTTCATTATTGTTTATTGGACATACTCTGT
lenght sequence : 96

sequence number :  9
start :  1485
frame :  1
sequence :  ATGGAAGGACAAAGTGCTGCGAAGTCATACGCTTCCAATTTCATTATTGTTTATTGGACATACTCTGTTAGCTTTATTACCGTCCACGCTTTTTCTACAATAGTG
lenght sequence : 105

sequence number :  10
start :  1830
frame :  1
sequence :  ATGGCTGGCTTTAATCTGCTGGAGTACCATGGAACACCGGTGATCATTCTGGTCACTTGGTCTGGGGCAATACCAGTCAACATGGTGGTGAAGTCACCG
lenght sequence : 99

sequence number :  11
start :  2107
frame :  2
sequence :  ATGGCAGCGACACCAGCAGCGATTGAAGTTAATTTGACCATTGTATTTGTTTTGTTTGTTAGTGCTGATGTAATCTTAACAAGAAATAGTGAAATGAAAGCGCATACCTCAAAGGCATATAGT
lenght sequence : 123

sequence number :  12
start :  2356
frame :  2
sequence :  ATGGAGATGCTCTTGTTTCTGAACGAATCATACATCTTTCATAGGTTTCGTATGTGGAGTATTGTTTTATGGCACTCATGTGTATTCGTATGCGCAGAATGTGGGAATGCCTATTATAGGGGTGCCGGGGGGTGCCTTGAAAAACCCTTTTGCGCGCCTGTTAAGTTTCCGTTTTCAGTCAAAAAGAATATCCGAATTTTAGATTTGGACCCTCGTTCAGAAGCTTACTGTCTAAGCCATCATTTGGTGTGTCCTAAACGGTTTCCATGCAAAGCAACATCTTTGTTACTCATTCCTGAAGGT
lenght sequence : 303

sequence number :  13
start :  2362
frame :  2
sequence :  ATGCTCTTGTTTCTGAACGAATCATACATCTTTCATAGGTTTCGTATGTGGAGTATTGTTTTATGGCACTCATGTGTATTCGTATGCGCAGAATGTGGGAATGCCTATTATAGGGGTGCCGGGGGGTGCCTTGAAAAACCCTTTTGCGCGCCTGTTAAGTTTCCGTTTTCAGTCAAAAAGAATATCCGAATTTTAGATTTGGACCCTCGTTCAGAAGCTTACTGTCTAAGCCATCATTTGGTGTGTCCTAAACGGTTTCCATGCAAAGCAACATCTTTGTTACTCATTCCTGAAGGT
lenght sequence : 297

sequence number :  14
start :  2407
frame :  2
sequence :  ATGTGGAGTATTGTTTTATGGCACTCATGTGTATTCGTATGCGCAGAATGTGGGAATGCCTATTATAGGGGTGCCGGGGGGTGCCTTGAAAAACCCTTTTGCGCGCCTGTTAAGTTTCCGTTTTCAGTCAAAAAGAATATCCGAATTTTAGATTTGGACCCTCGTTCAGAAGCTTACTGTCTAAGCCATCATTTGGTGTGTCCTAAACGGTTTCCATGCAAAGCAACATCTTTGTTACTCATTCCTGAAGGT
lenght sequence : 252

sequence number :  15
start :  2462
frame :  3
sequence :  ATGCCTATTATAGGGGTGCCGGGGGGTGCCTTGAAAAACCCTTTTGCGCGCCTGTTAAGTTTCCGTTTTCAGTCAAAAAGAATATCCGAATTT
lenght sequence : 93

sequence number :  16
start :  3447
frame :  1
sequence :  ATGATAATAATATTTATAGAATTGTGTAGAATTGCAGATTCCCTTTTATGGATTCCTAAATCCTGGAGGAGAACTTCTAGTATATCTACATACCTAATATTA
lenght sequence : 102

sequence number :  17
start :  3563
frame :  3
sequence :  ATGGAATCCCAACAATTACATCAAAATCCACATTCTCTACACCAATACCATCGACGAGAGCTTCTAGTAAATTGTATACATAACAGTATAACCCTTACCAACAATGGAATCTCAAAGATTATTAAATTATTCACAGACTCTGAGGATTCGGGTAAAATAGGGTATTTAACTGGTTACCGGAAAGGTTTAGAAAATTCGTGGAGGGTTGGCCGAGTGGTC
lenght sequence : 219

sequence number :  18
start :  4056
frame :  1
sequence :  ATGATGTGTAGCGCATTAGTCTTCTGCATTCATATTATCAGTGGTATACCGTTTGGTATATGCTTTCCAATCATGAAGTGTTGGCTCTACTCTATA
lenght sequence : 96

sequence number :  19
start :  4059
frame :  1
sequence :  ATGTGTAGCGCATTAGTCTTCTGCATTCATATTATCAGTGGTATACCGTTTGGTATATGCTTTCCAATCATGAAGTGTTGGCTCTACTCTATA
lenght sequence : 93

sequence number :  20
start :  5073
frame :  1
sequence :  ATGCTGCGACAGAACTTTTCTTCTCTTCCAATTCTTTTATCTTCGCCAATAACGCCTTCTTTTTTCGTGCTTGTATTTCTAAAATTTCGGCCAGGTATTGTAAATATTCGACCGTTCTATCCAAAA
lenght sequence : 126

sequence number :  21
start :  5253
frame :  1
sequence :  ATGGTGGAACCAACTGGCCGAGTTCTTTTATCTTTTGCTTTATTAATTCTCTTCTTCTCCTTTCGACGGCATTATGAAACTCTCTTTTGCGCCTCAGTTTCTCATCAGAAGTTAACCCGCCTAAAATCTTTGGAACACTTCCAGGTCCTATATTTTCAGTCATATTGCTACTTATTGAAGTGTGTCTTGTTCTGGGTGTGTTTATGCTACCATGCC
lenght sequence : 216

sequence number :  22
start :  5456
frame :  3
sequence :  ATGCTACCATGCCTAAAAGAAGATGAAAGGAAACTTCCTGCACGAAATGACGATGATGGTGATCTTACCTTTGGGGAATATGTGGAAGATACGGATGCTGGGCCCAAGCTTTGTGGATTA
lenght sequence : 120

sequence number :  23
start :  5478
frame :  1
sequence :  ATGAAAGGAAACTTCCTGCACGAAATGACGATGATGGTGATCTTACCTTTGGGGAATATGTGGAAGATACGGATGCTGGGCCCAAGCTTTGTGGATTATAAGAAAATGATGATGAATATGTGTTTGGTGTCATCATATCAGAAT
lenght sequence : 144

sequence number :  24
start :  5502
frame :  1
sequence :  ATGACGATGATGGTGATCTTACCTTTGGGGAATATGTGGAAGATACGGATGCTGGGCCCAAGCTTTGTGGATTATAAGAAAATGATGATGAATATGTGTTTGGTGTCATCATATCAGAAT
lenght sequence : 120

sequence number :  25
start :  5508
frame :  1
sequence :  ATGATGGTGATCTTACCTTTGGGGAATATGTGGAAGATACGGATGCTGGGCCCAAGCTTTGTGGATTATAAGAAAATGATGATGAATATGTGTTTGGTGTCATCATATCAGAAT
lenght sequence : 114

sequence number :  26
start :  5511
frame :  1
sequence :  ATGGTGATCTTACCTTTGGGGAATATGTGGAAGATACGGATGCTGGGCCCAAGCTTTGTGGATTATAAGAAAATGATGATGAATATGTGTTTGGTGTCATCATATCAGAAT
lenght sequence : 111

I am sure there is some awk, perl, regex way to get the same in much less lines.

Thank you for your edit

I am gonna post another answer using regex copied from https://stackoverflow.com/questions/13114246/how-to-find-a-open-reading-frame-in-python please feel free to edit it solving the start position problem

Thanks @pippo1980 good you've done this thanks upvoted again. The protein is now identified and a better approach IMO is to map against the genome.

using an unwanted (regex library) I can steal code from https://stackoverflow.com/questions/13114246/how-to-find-a-open-reading-frame-in-python , code :
import re

from Bio import SeqIO
# record = SeqIO.read("sample_test.fasta", "fasta")
# record = SeqIO.read("sample_test_2.fasta", "fasta")

record = SeqIO.read("sample.fasta", "fasta")

# print(record.seq)
# print(len(sequence))

sequence = str(record.seq)

pattern = re.compile(r'(?=(ATG(?:...)*?)(?=TAG|TGA|TAA))')

revcompseq = sequence[::-1].maketrans("ATGC", "TACG") #reverse complement

# print (pattern.findall(sequence)) #forward search
# print (pattern.findall(sequence[::-1].translate(revcompseq))) #backward search

a = pattern.findall(sequence)

# print(type(a))

cnt = 1
for i in a:
    if len(i) > 90:
        print('\nnumero = ',cnt,'\n',i ,'\nlenght = ', len(i),'\n---------------------------\n\n')
        cnt += 1

output :
numero =  1
 ATGCAATACAGTGACATATATATATACCCTAACACTACCCTAACCCTACCCTATTTCAACCCTTCCAACCTGTCTCTCAACTTACCCTCACATTACCCTACCTCTCCACTTGTTACCCTGTCCCATTCAACCATACCACTCCCAACCACCATCCATCCCTCTACTTACTACCACCAATCAACCGTCCACCATAACCGTTACCCTCCAATTAGCCATATTCAACTTCACTACCACTTACCCTGCCATTACTCTACCATCCACCATCTGCTACTCACCATACTGTTGTTCTACCCTCCATAT
lenght =  300
---------------------------

numero =  2
 ATGCACGGCACTTGCCTCAGCGGTCTATACCCTGAGCCATTTACCCATAACTCCCACGATTATCCACATTTTAATATCTATATCTCATTCGGCGGGCCCAAATATTGTATAACTGCTCTTAATACATACGTTATACCACTTTTGCACCATATACTAACCACTCAAATTATATACACTTATGCCAATATAACCAAAAAATTACCACTAAAATCACCTAAACATAAAAATTATTTATCTTTCAACTTTACGAAA
lenght =  252
---------------------------

numero =  3
 ATGCCACCTCACTGTCGTAACACTCTATATTCACCGAGAAACAATACGGTAGTGGCTCAAACTCATGCCGGTGCTATGATACAATTGTATCCTATTTCCATTCTCATATGCTATCCGCAATATCCTAAAAGCATAACTGATGCATCTTTAATCTTGTATGTGACACTACTCATACGAAGAGACTATATC
lenght =  189
---------------------------

numero =  4
 ATGATACAATTGTATCCTATTTCCATTCTCATATGCTATCCGCAATATCCTAAAAGCATAACTGATGCATCTTTAATCTTGTATGTGACACTACTCATACGAAGAGACTATATC
lenght =  114
---------------------------

numero =  5
 ATGAAACGTCAAATAACCCTACTATATTATAACTTATCAGCGCCGTATACTAAAACGGACGTTACGATATTGTCTCACTTCATCTTACCACCCTCTATCTTATTGTTGATAAAACAC
lenght =  117
---------------------------

numero =  6
 ATGACTTGTAACTCCACTACCCTGATCTGCAATCTTGTTCTTAGAAGTGACGCATACTCTATATGGCCCGACGCGACGCGCCAAAAAATGAAAAAAGAAGCAGCGACTCATTTTTATGGAAGGACAAAGTGCTGCGAAGTCATACGCTTCCAATTTCATTATTGTTTATTGGACATACTCTGT
lenght =  183
---------------------------

numero =  7
 ATGGCCCGACGCGACGCGCCAAAAAATGAAAAAAGAAGCAGCGACTCATTTTTATGGAAGGACAAAGTGCTGCGAAGTCATACGCTTCCAATTTCATTATTGTTTATTGGACATACTCTGTTAGCTTTATTACCGTCCACGCTTTTTCTACAA
lenght =  153
---------------------------

numero =  8
 ATGAAAAAAGAAGCAGCGACTCATTTTTATGGAAGGACAAAGTGCTGCGAAGTCATACGCTTCCAATTTCATTATTGTTTATTGGACATACTCTGT
lenght =  96
---------------------------

numero =  9
 ATGGAAGGACAAAGTGCTGCGAAGTCATACGCTTCCAATTTCATTATTGTTTATTGGACATACTCTGTTAGCTTTATTACCGTCCACGCTTTTTCTACAATAGTG
lenght =  105
---------------------------

numero =  10
 ATGGCTGGCTTTAATCTGCTGGAGTACCATGGAACACCGGTGATCATTCTGGTCACTTGGTCTGGGGCAATACCAGTCAACATGGTGGTGAAGTCACCG
lenght =  99
---------------------------

numero =  11
 ATGGCAGCGACACCAGCAGCGATTGAAGTTAATTTGACCATTGTATTTGTTTTGTTTGTTAGTGCTGATGTAATCTTAACAAGAAATAGTGAAATGAAAGCGCATACCTCAAAGGCATATAGT
lenght =  123
---------------------------

numero =  12
 ATGGAGATGCTCTTGTTTCTGAACGAATCATACATCTTTCATAGGTTTCGTATGTGGAGTATTGTTTTATGGCACTCATGTGTATTCGTATGCGCAGAATGTGGGAATGCCTATTATAGGGGTGCCGGGGGGTGCCTTGAAAAACCCTTTTGCGCGCCTGTTAAGTTTCCGTTTTCAGTCAAAAAGAATATCCGAATTTTAGATTTGGACCCTCGTTCAGAAGCTTACTGTCTAAGCCATCATTTGGTGTGTCCTAAACGGTTTCCATGCAAAGCAACATCTTTGTTACTCATTCCTGAAGGT
lenght =  303
---------------------------

numero =  13
 ATGCTCTTGTTTCTGAACGAATCATACATCTTTCATAGGTTTCGTATGTGGAGTATTGTTTTATGGCACTCATGTGTATTCGTATGCGCAGAATGTGGGAATGCCTATTATAGGGGTGCCGGGGGGTGCCTTGAAAAACCCTTTTGCGCGCCTGTTAAGTTTCCGTTTTCAGTCAAAAAGAATATCCGAATTTTAGATTTGGACCCTCGTTCAGAAGCTTACTGTCTAAGCCATCATTTGGTGTGTCCTAAACGGTTTCCATGCAAAGCAACATCTTTGTTACTCATTCCTGAAGGT
lenght =  297
---------------------------

numero =  14
 ATGTGGAGTATTGTTTTATGGCACTCATGTGTATTCGTATGCGCAGAATGTGGGAATGCCTATTATAGGGGTGCCGGGGGGTGCCTTGAAAAACCCTTTTGCGCGCCTGTTAAGTTTCCGTTTTCAGTCAAAAAGAATATCCGAATTTTAGATTTGGACCCTCGTTCAGAAGCTTACTGTCTAAGCCATCATTTGGTGTGTCCTAAACGGTTTCCATGCAAAGCAACATCTTTGTTACTCATTCCTGAAGGT
lenght =  252
---------------------------

numero =  15
 ATGCCTATTATAGGGGTGCCGGGGGGTGCCTTGAAAAACCCTTTTGCGCGCCTGTTAAGTTTCCGTTTTCAGTCAAAAAGAATATCCGAATTT
lenght =  93
---------------------------

numero =  16
 ATGATAATAATATTTATAGAATTGTGTAGAATTGCAGATTCCCTTTTATGGATTCCTAAATCCTGGAGGAGAACTTCTAGTATATCTACATACCTAATATTA
lenght =  102
---------------------------

numero =  17
 ATGGAATCCCAACAATTACATCAAAATCCACATTCTCTACACCAATACCATCGACGAGAGCTTCTAGTAAATTGTATACATAACAGTATAACCCTTACCAACAATGGAATCTCAAAGATTATTAAATTATTCACAGACTCTGAGGATTCGGGTAAAATAGGGTATTTAACTGGTTACCGGAAAGGTTTAGAAAATTCGTGGAGGGTTGGCCGAGTGGTC
lenght =  219
---------------------------

numero =  18
 ATGATGTGTAGCGCATTAGTCTTCTGCATTCATATTATCAGTGGTATACCGTTTGGTATATGCTTTCCAATCATGAAGTGTTGGCTCTACTCTATA
lenght =  96
---------------------------

numero =  19
 ATGTGTAGCGCATTAGTCTTCTGCATTCATATTATCAGTGGTATACCGTTTGGTATATGCTTTCCAATCATGAAGTGTTGGCTCTACTCTATA
lenght =  93
---------------------------

numero =  20
 ATGCTGCGACAGAACTTTTCTTCTCTTCCAATTCTTTTATCTTCGCCAATAACGCCTTCTTTTTTCGTGCTTGTATTTCTAAAATTTCGGCCAGGTATTGTAAATATTCGACCGTTCTATCCAAAA
lenght =  126
---------------------------

numero =  21
 ATGGTGGAACCAACTGGCCGAGTTCTTTTATCTTTTGCTTTATTAATTCTCTTCTTCTCCTTTCGACGGCATTATGAAACTCTCTTTTGCGCCTCAGTTTCTCATCAGAAGTTAACCCGCCTAAAATCTTTGGAACACTTCCAGGTCCTATATTTTCAGTCATATTGCTACTTATTGAAGTGTGTCTTGTTCTGGGTGTGTTTATGCTACCATGCC
lenght =  216
---------------------------

numero =  22
 ATGCTACCATGCCTAAAAGAAGATGAAAGGAAACTTCCTGCACGAAATGACGATGATGGTGATCTTACCTTTGGGGAATATGTGGAAGATACGGATGCTGGGCCCAAGCTTTGTGGATTA
lenght =  120
---------------------------

numero =  23
 ATGAAAGGAAACTTCCTGCACGAAATGACGATGATGGTGATCTTACCTTTGGGGAATATGTGGAAGATACGGATGCTGGGCCCAAGCTTTGTGGATTATAAGAAAATGATGATGAATATGTGTTTGGTGTCATCATATCAGAAT
lenght =  144
---------------------------

numero =  24
 ATGACGATGATGGTGATCTTACCTTTGGGGAATATGTGGAAGATACGGATGCTGGGCCCAAGCTTTGTGGATTATAAGAAAATGATGATGAATATGTGTTTGGTGTCATCATATCAGAAT
lenght =  120
---------------------------

numero =  25
 ATGATGGTGATCTTACCTTTGGGGAATATGTGGAAGATACGGATGCTGGGCCCAAGCTTTGTGGATTATAAGAAAATGATGATGAATATGTGTTTGGTGTCATCATATCAGAAT
lenght =  114
---------------------------

numero =  26
 ATGGTGATCTTACCTTTGGGGAATATGTGGAAGATACGGATGCTGGGCCCAAGCTTTGTGGATTATAAGAAAATGATGATGAATATGTGTTTGGTGTCATCATATCAGAAT
lenght =  111
---------------------------

I can get the start position in my sequence of all the found ORFs with :
b = [[m.start(1), m.end(1)] for m in pattern.finditer(sequence)]

print(b, len(b))

result a list of tuples :
[(2, 302), (316, 400), (373, 400), (404, 479), (421, 505),
(478, 505), (508, 556), (524, 776), (702, 711), (807, 822),
(907, 1096), (971, 983), (982, 1096), (1014, 1032), (1046, 1055),
(1064, 1067), (1117, 1141), (1143, 1260), (1340, 1355), (1354, 1378),
(1370, 1553), (1432, 1585), (1457, 1553), (1485, 1590), (1606, 1663),
 (1626, 1647), (1754, 1769), (1830, 1929),

.........................................

but cannot filter for m.end() - m.start() > 90 in the one liner b
any help ??
EDITED :
b = [[m.start(1), m.end(1)] for m in pattern.finditer(sequence) if (m.end(1) - m.start(1)) > 90]

should work only cant get found string with m.group()
as an example for first match I get:
2, 302, <re.Match object; span=(2, 2), match=''> with m.group() = ''
EDITED 2:
I can get the found substring using  m.groups()[0]

The code is,
from Bio.Seq import Seq
frame = 2
myseq = myseq[frame:]
bioseq = Seq(myseq)
fullaa = bioseq.translate()
for aa in fullaa.split("*"):
    print (aa)

The sequence is below the -----
The output looks like open readingframe2. The first protein is:

MQYSDIYIYPNTTLTLPYFNPSNLSLNLPSHYPTSPLVTLSHSTIPLPTTIHPSTYYHQSTVHHNRYPPISHIQLHYHLPCHYSTIHHLLLTILLFYPPY

This meets with the requirements ... the remaining output is below ... it looks like a eukaryotic with introns. It feels like a home-work question.

Edit Blasted it, its a yeast protein

RecName: Full=Putative uncharacterized protein YBL109W [Saccharomyces cerevisiae S288C]
Sequence ID: P89495.3

The hit is 100%

Remaining output
NVNK
S
IIHIYLPYHFIPPHITCHTHLHLYTDMPYAHGCYSIYHSQTYPTLTFYSTP
PISH
ISTKCTHIIMHGTCLSGLYPEPFTHNSHDYPHFNIYISFGGPKYCITALNTYVIPLLHHILTTQIIYTYANITKKLPLKSPKHKNYLSFNFTK
THSIAYLYHHDVINVKVP
YYHLLERIPFQNISNFHRPYIRIICHLTVVTLYIHRETIR
WLKLMPVL
YNCILFPFSYAIRNILKA
LMHL
SCM
HYSYEETISKEDDTVICTLFL
NYNETSNNPTIL
LISAVY
NGRYDIVSLHLTTLYLIVDKTLTPQLYF
LQLQQTIPNHKS
YFRCQKMRISK
EFGTMTCNSTTLICNLVLRSDAYSIWPDATRQKMKKEAATHFYGRTKCCEVIRFQFHYCLLDILC
LYYRPRFFYNSVKVSFLCSSYS
NASRTPSLIK
VYNINIHLYNRRYLYHQKK
FFILFFHYFSLSMDFHS
RHHSLVCDSVDTVLG
STGDGWL
SAGVPWNTGDHSGHLVWGNTSQHGGEVTVVENGFSNFDWVGFSWVGGLEHVVLGQVSSDIRDVDTQFHQVDSFVRLS
SGGCRKQ
QRWQRHQQRLKLI
PLYLFCLLVLM
S
QEIVK
KRIPQRHIVEAALFIPVPPSVITT
TIR
QTLI
HLTYSPYLIFHTISLSLWRCSCF
TNHTSFIGFVCGVLFYGTHVYSYAQNVGMPIIGVPGGALKNPFARLLSFRFQSKRISEF
IWTLVQKLTV
AIIWCVLNGFHAKQHLCYSFLKVE
KIATMCRESYIVRF
RLVRIIKAVNRINIRKASKVHS
ENS
HVKICVYEMCQGPSIA
LTPLQE
VTFFTLKRQRGHILKSHFQTHLQFATTDNNIMIGWYDVISCHIFLIIALVTSLQLLG
NSKISFL

YSSYQYIITCRKRR
HKDRQTVFKYNGSWNARIDNVTGY
MTKYK
KKKKKKSSNTIMWKYRFHFEDSNCWNKNQLSSTN
YLRY
YIIIYGVRR
RK

EIVI
ISGS
NARIDNVIGSMNINI
NDDNNIYRIV
NCRFPFMDS
ILEENF
YIYIPNIIALITMESQQLHQNPHSLHQYHRRELLVNCIHNSITLTNNGISKIIKLFTDSEDSGKIGYLTGYRKGLENSWRVGRVV
GGRLKICWTVVRASSNLAFIFRIHFIVGHFISCRVH
KIINQIII
PMGHFFDETHVYSSIRS
YCTYDYQ
HQNNHRAYAILKYIIAVKYDV
RISLLHSYYQWYTVWYMLSNHEVLALLYINSTSIDFLNK
VTQVSFPLIQVIYILATPNYKALKIISYPIAGSTCIYYTLLITNLLHI
QYA

Y
CI
IIS
KKITWVLLT
LHLEYIKF
RINISKVLDITKNNSE
F
ILLHIILLAYLFHSFQNVACKRY
YSALFCCNS
FFHYIILRVRKRLTFSQ
KEQT
ILKFYCFCLPDLFQI
FFPANKTINYPEPNSKR
SAGLELSSPNINFALSSASIKSPDKNSFNSS
MFVCG
LSIVPA
LWLPTESYNGGSQLCNSALLFCSFIKAFGTDLISIILSSLFSLLFCPELAWWLLVKGDNAATELFFSSNSFIFANNAFFFRACISKISARYCKYSTVLSKMIPLFGLICLPRSS
FNKDGGTNWPSSFIFCFINSLLLLSTAL
NSLLRLSFSSEVNPPKIFGTLPGPIFSVILLLIEVCLVLGVFMLPCLKEDERKLPARNDDDGDLTFGEYVEDTDAGPKLCGL
ENDDEYVFGVIISELMLEDKEELKSSVKLYIVSSVPC
NKKSCGNSAFKPSCCCNGVFIA

myseq = "GCATGCAATACAGTGACATATATATATACCCTAACACTACCCTAACCCTACCCTATTTCA\
ACCCTTCCAACCTGTCTCTCAACTTACCCTCACATTACCCTACCTCTCCACTTGTTACCC\
TGTCCCATTCAACCATACCACTCCCAACCACCATCCATCCCTCTACTTACTACCACCAAT\
CAACCGTCCACCATAACCGTTACCCTCCAATTAGCCATATTCAACTTCACTACCACTTAC\
CCTGCCATTACTCTACCATCCACCATCTGCTACTCACCATACTGTTGTTCTACCCTCCAT\
ATTGAAACGTTAACAAATGATCGTAAATAATACACATATACTTACCCTACCACTTTATAC\
CACCACACATCACATGCCATACTCACCTTCACTTGTATACTGATATGCCATACGCACACG\
GATGCTACAGTATATACCACTCTCAAACTTACCCTACTCTCACATTCTACTCCACTCCAT\
GACCCATCTCTCACTAAATCAGTACTAAATGCACCCACATCATTATGCACGGCACTTGCC\
TCAGCGGTCTATACCCTGAGCCATTTACCCATAACTCCCACGATTATCCACATTTTAATA\
TCTATATCTCATTCGGCGGGCCCAAATATTGTATAACTGCTCTTAATACATACGTTATAC\
CACTTTTGCACCATATACTAACCACTCAAATTATATACACTTATGCCAATATAACCAAAA\
AATTACCACTAAAATCACCTAAACATAAAAATTATTTATCTTTCAACTTTACGAAATAAA\
CACACTCAATTGCGTATCTATACCACCATGACGTCATTAACGTAAAAGTTCCTTAATATT\
ACCATTTGCTTGAACGGATACCATTTCAGAATATTTCTAACTTTCACAGACCATACATTA\
GAATAATATGCCACCTCACTGTCGTAACACTCTATATTCACCGAGAAACAATACGGTAGT\
GGCTCAAACTCATGCCGGTGCTATGATACAATTGTATCCTATTTCCATTCTCATATGCTA\
TCCGCAATATCCTAAAAGCATAACTGATGCATCTTTAATCTTGTATGTGACACTACTCAT\
ACGAAGAGACTATATCTAAAGAAGACGATACAGTGATATGTACGTTGTTTTTGTAGAATT\
ATAATGAAACGTCAAATAACCCTACTATATTATAACTTATCAGCGCCGTATACTAAAACG\
GACGTTACGATATTGTCTCACTTCATCTTACCACCCTCTATCTTATTGTTGATAAAACAC\
TAACCCCTCAGCTTTATTTCTAGTTACAGTTACAACAAACTATCCCAAACCATAAATCTT\
AATATTTTAGGTGTCAAAAAATGAGGATCTCCAAATGAGAGTTTGGTACCATGACTTGTA\
ACTCCACTACCCTGATCTGCAATCTTGTTCTTAGAAGTGACGCATACTCTATATGGCCCG\
ACGCGACGCGCCAAAAAATGAAAAAAGAAGCAGCGACTCATTTTTATGGAAGGACAAAGT\
GCTGCGAAGTCATACGCTTCCAATTTCATTATTGTTTATTGGACATACTCTGTTAGCTTT\
ATTACCGTCCACGCTTTTTCTACAATAGTGTAAAAGTTTCTTTCTTATGTTCATCGTATT\
CATAAAATGCTTCACGAACACCGTCATTGATCAAATAGGTTTATAATATTAATATACATT\
TATATAATCGGCGGTATTTATATCATCAAAAAAAGTAGTTTTTTATTTTATTTTTTCATT\
ACTTTTCACTGTCTATGGATTTTCATTCGTAAAGGCATCACTCCCTAGTTTGCGATAGTG\
TAGATACCGTCCTTGGATAGAGCACTGGAGATGGCTGGCTTTAATCTGCTGGAGTACCAT\
GGAACACCGGTGATCATTCTGGTCACTTGGTCTGGGGCAATACCAGTCAACATGGTGGTG\
AAGTCACCGTAGTTGAAAACGGCTTCAGCAACTTCGACTGGGTAGGTTTCAGTTGGGTGG\
GCGGCTTGGAACATGTAGTATTGGGCCAAGTGAGCTCTGATATCAGAGACGTAGACACCC\
AATTCCACCAAGTTGACTCTTTCGTCAGATTGAGCTAGAGTGGTGGTTGCAGAAAGCAGT\
AGCAGCGATGGCAGCGACACCAGCAGCGATTGAAGTTAATTTGACCATTGTATTTGTTTT\
GTTTGTTAGTGCTGATGTAATCTTAACAAGAAATAGTGAAATGAAAGCGCATACCTCAAA\
GGCATATAGTTGAAGCAGCTCTATTTATACCCGTTCCTCCATCTGTCATCACTACTTAAA\
CGATTCGTTAACAGACGCTCATTTAGCACCTCACATATTCTCCATATCTCATCTTTCACA\
CAATCTCATTATCTCTATGGAGATGCTCTTGTTTCTGAACGAATCATACATCTTTCATAG\
GTTTCGTATGTGGAGTATTGTTTTATGGCACTCATGTGTATTCGTATGCGCAGAATGTGG\
GAATGCCTATTATAGGGGTGCCGGGGGGTGCCTTGAAAAACCCTTTTGCGCGCCTGTTAA\
GTTTCCGTTTTCAGTCAAAAAGAATATCCGAATTTTAGATTTGGACCCTCGTTCAGAAGC\
TTACTGTCTAAGCCATCATTTGGTGTGTCCTAAACGGTTTCCATGCAAAGCAACATCTTT\
GTTACTCATTCCTGAAGGTTGAATAAAAAATAGCAACCATGTGCAGGGAGTCGTATATTG\
TTAGATTCTAGAGACTTGTACGCATCATCAAAGCTGTAAATAGAATAAACATACGCAAGG\
CGTCAAAAGTGCATAGTTAAGAAAATTCCTGACATGTGAAAATATGTGTTTATGAAATGT\
GTCAAGGCCCGTCTATAGCGTAGTTAACCCCTCTGCAGGAGTAAGTGACTTTTTTTACGC\
TCAAAAGGCAACGAGGGCACATACTTAAAAGTCATTTTCAAACACATCTGCAGTTTGCAA\
CGACAGATAACAATATTATGATAGGATGGTATGATGTTATTAGCTGCCACATATTTTTAA\
TAATAGCATTAGTCACGTCTCTTCAATTGTTGGGATGAAACTCTAAAATATCATTCCTTT\
AGTAGTATTCCAGTTACCAGTATATTATCACATGCCGAAAAAGAAGATGACATAAAGATC\
GACAAACAGTCTTCAAATATAATGGAAGCTGGAATGCAAGGATTGATAATGTAACAGGAT\
ACTGAATGACAAAGTATAAATGAAAAAAAAAAAAAAAAAAAAGTAGTAATACTATTATGT\
GGAAATACCGATTCCATTTTGAGGATTCCAATTGTTGGAATAAAAATCAACTATCATCTA\
CTAACTAGTATTTACGTTACTAGTATATTATCATATACGGTGTTAGAAGATGACGCAAAT\
GATGAGAAATAGTCATCTAAATTAGTGGAAGCTGAAACGCAAGGATTGATAATGTAATAG\
GATCAATGAATATTAACATATAAAACGATGATAATAATATTTATAGAATTGTGTAGAATT\
GCAGATTCCCTTTTATGGATTCCTAAATCCTGGAGGAGAACTTCTAGTATATCTACATAC\
CTAATATTATAGCCTTAATCACAATGGAATCCCAACAATTACATCAAAATCCACATTCTC\
TACACCAATACCATCGACGAGAGCTTCTAGTAAATTGTATACATAACAGTATAACCCTTA\
CCAACAATGGAATCTCAAAGATTATTAAATTATTCACAGACTCTGAGGATTCGGGTAAAA\
TAGGGTATTTAACTGGTTACCGGAAAGGTTTAGAAAATTCGTGGAGGGTTGGCCGAGTGG\
TCTAAGGCGGCAGACTTAAGATCTGTTGGACGGTTGTCCGCGCGAGTTCGAACCTCGCAT\
TCATTTTCCGGATCCACTTTATAGTGGGTCACTTCATCTCTTGTAGGGTTCACTGAAAGA\
TAATCAATCAAATTATCATATGACCAATGGGTCACTTTTTTGATGAGACCCATGTTTATT\
CTTCTATACGTTCGTGATACTGCACTTACGACTACCAGTAACATCAGAATAACCATAGAG\
CTTATGCTATCTTGAAATATATCATCGCCGTAAAGTATGATGTGTAGCGCATTAGTCTTC\
TGCATTCATATTATCAGTGGTATACCGTTTGGTATATGCTTTCCAATCATGAAGTGTTGG\
CTCTACTCTATATAAATTCGACATCCATTGATTTTCTCAATAAATGAGTTACCCAAGTAA\
GCTTTCCATTGATACAAGTGATCTACATTCTTGCGACGCCAAATTATAAAGCACTAAAAA\
TCATATCATACCCAATTGCGGGCAGCACATGTATATATTATACACTACTTATAACTAACC\
TTCTTCATATATAACAATATGCCTGATAATACTAATGCATTTAAATCATATCATGAAAGA\
AAATTACATGGGTTTTATTGACATAATTGCATTTAGAATACATAAAATTCTAAAGAATTA\
ATATATCCAAAGTATTAGACATAACCAAGAATAATAGTGAATAATTTTAGATTTTGTTAC\
ATATAATTCTGCTTGCCTATCTCTTCCACTCTTTTCAAAACGTTGCATGTAAGCGTTACT\
AATATTCCGCTTTATTTTGTTGCAATTCCTAATTTTTTCATTACATTATCTTGCGAGTAC\
GGAAGCGATTAACGTTCTCCCAATAGAAGGAACAAACATAGATATTGAAGTTTTACTGCT\
TTTGCTTACCTGACCTTTTTCAAATTTAATTTTTTCCCGCTAATAAGACCATAAACTACC\
CCGAACCAAATTCTAAAAGATAGTCAGCTGGATTAGAGTTGTCATCTCCAAACATTAATT\
TTGCATTATCTTCGGCTTCAATCAAATCGCCTGATAAGAACTCCTTTAATTCTTCATGAA\
TGTTTGTATGTGGATGACTCTCCATAGTGCCAGCATGATTGTGGTTACCGACCGAATCAT\
ATAACGGTGGCTCCCAATTGTGCAATTCAGCCTTACTATTTTGTTCATTCATCAAAGCAT\
TTGGGACAGATCTAATATCTATAATTCTTTCCTCACTATTCTCGCTATTATTTTGCCCCG\
AACTGGCATGGTGGTTATTGGTAAAAGGAGATAATGCTGCGACAGAACTTTTCTTCTCTT\
CCAATTCTTTTATCTTCGCCAATAACGCCTTCTTTTTTCGTGCTTGTATTTCTAAAATTT\
CGGCCAGGTATTGTAAATATTCGACCGTTCTATCCAAAATGATACCCTTATTTGGTTTGA\
TTTGTTTACCTAGGTCATCGTAATTCAATAAAGATGGTGGAACCAACTGGCCGAGTTCTT\
TTATCTTTTGCTTTATTAATTCTCTTCTTCTCCTTTCGACGGCATTATGAAACTCTCTTT\
TGCGCCTCAGTTTCTCATCAGAAGTTAACCCGCCTAAAATCTTTGGAACACTTCCAGGTC\
CTATATTTTCAGTCATATTGCTACTTATTGAAGTGTGTCTTGTTCTGGGTGTGTTTATGC\
TACCATGCCTAAAAGAAGATGAAAGGAAACTTCCTGCACGAAATGACGATGATGGTGATC\
TTACCTTTGGGGAATATGTGGAAGATACGGATGCTGGGCCCAAGCTTTGTGGATTATAAG\
AAAATGATGATGAATATGTGTTTGGTGTCATCATATCAGAATTGATGCTAGAAGATAAAG\
AGGAGCTCAAATCATCAGTTAAGTTATACATTGTGTCATCGGTGCCATGCTGAAATAAAA\
AGTCATGTGGTAACTCAGCTTTTAAGCCTTCTTGTTGCTGTAATGGAGTATTCATTGCTC"


text:
Creating meta cells from scRNA-seq and scATAC-seq for GRN analysis
We have two 10X Genomics scRNA-seq and scATAC-seq datasets generated from similar cells, but in different experiments - separately, not using the multiome platform. Now, we are trying to use Pando (https://github.com/quadbiolab/Pando) to generate a gene regulatory network.
We have successfully converted our CCA-embedded dataset into h5ad format (https://mojaveazure.github.io/seurat-disk/articles/convert-anndata.html), then used bipartite matching (https://github.com/quadbiolab/organoid_regulomes) to receive a .tsv file matching invidiual cells from both datasets - it gives us the barcodes of the cells in both the scRNA and scATAC datasets. Now, we need to actually generate metacells/pseudocells that contain both modalities inside the same cell to be able to use Pando.
According to Fleck et al. 2022 (https://www.nature.com/articles/s41586-022-05279-8), they state

The metacells were integrated using cluster similarity spectrum (CSS)

which would refer to this technique: https://github.com/quadbiolab/simspec/blob/master/vignette/vignette.md
However, that does not make use of the previously generated bipartite-matched cell list.
I'm a little bit lost here - can someone point me in the right direction to how we can generate the metacells containing both modalities, scRNA and scATAC?
Thanks a lot in advance!

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
Trim reads 1kb upstream of sequence
I need a quick way to trim multiple reads in a FASTA file. I need to trim everything that is 1kbp upstream of this sequence AAGAGATGTTCAATCGTTTAAACAAATTCCAAGCTGCTTTAGCTTTGGCCCTTTACTCTCA.
I figure a quick python script might be the way to go but I'm not sure if there's a tool that already does that (i.e. Trimmomatic or something like it) that would be able to do that for me.
Thank you!

What is the source of these reads (i.e. sequencing platform)? Is it necessary to account for some error in the recognition sequence?

Here's one way using biopython. Note that the Seq object has a number of methods that act just like those of a Python string. One of these is the find() method, which returns the index of the first occurrence of the sequence if found. Replace with rfind() if you want the last occurrence if that makes more sense. This solution uses a generator to avoid storing the entire list of trimmed sequences in memory:
import argparse
import contextlib
import gzip
import sys

from Bio import SeqIO

def trim_fasta(fileobj, sequence, num_bases):
    """
    Trim num_bases upstream of the sequence provided
    """
    for record in SeqIO.parse(fileobj, 'fasta'):
        index = record.seq.find(sequence)
        if index != -1 and index > num_bases:
            record = record[index-num_bases:]
        yield record

def open_file(filename, mode='r'):
    if filename.endswith('.gz'):
        return gzip.open(filename, mode)
    else:
        return open(filename, mode)

def get_argument_parser():

    parser = argparse.ArgumentParser(add_help=False)

    group = parser.add_argument_group('input options')
    group.add_argument(
        "fasta",
        type=str,
        metavar="FILE",
        help="The input FASTA file",
    )
    group.add_argument(
        "sequence",
        type=str,
        metavar="STR",
        help="The sequence to find and trim",
    )

    group = parser.add_argument_group('trimming options')
    group.add_argument(
        "-n",
        "--num_bases",
        type=int,
        metavar="INT",
        default=1000,
        help="Trim INT bases upstream of input sequence",
    )

    group = parser.add_argument_group('output options')
    group.add_argument(
        "-o",
        "--output",
        type=str,
        metavar="FILE",
        default='-',
        help="Write the filtered output to FILE (default: stdout)",
    )
    group.add_argument(
        "-f",
        "--force",
        action='store_true',
        help="Overwrite the output file if it exists",
    )

    group = parser.add_argument_group('additional options')
    group.add_argument(
        "-h",
        "--help",
        action="help",
        help="Show this help message and exit",
    )

    return parser

def main():

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

    with contextlib.ExitStack() as stack:
        if args.output != '-':
            out = stack.enter_context(open_file(args.output, 'wt' if args.force else 'xt'))
        else:
            out = sys.stdout

        input_fasta = stack.enter_context(open_file(args.fasta))
        trimmed_reads = trim_fasta(input_fasta, args.sequence, args.num_bases)

        SeqIO.write(trimmed_reads, out, "fasta-2line")

if __name__ == '__main__':
    main()

usage: trim_reads.py [-n INT] [-o FILE] [-f] [-h] FILE STR

input options:
  FILE                  The input FASTA file
  STR                   The sequence to find and trim

trimming options:
  -n INT, --num_bases INT
                        Trim INT bases upstream of input sequence

output options:
  -o FILE, --output FILE
                        Write the filtered output to FILE (default: stdout)
  -f, --force           Overwrite the output file if it exists

additional options:
  -h, --help            Show this help message and exit

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
PCA of methylation-values normalized by coverage
I have started analyzing methylation (EM-seq) data for the first time ~0.8M positions.
In this data set I have 27 samples of 20 patients. I want to perform a PCA of the dataset to check for possible batch effects or problems.
Some publications and guides use the raw M-values to plot it without taking into account the coverage of each position. If I do that, I can detect a previously unknown source of variation, which is already great.
Usually there is also the information of the coverage of each position (not the number of reads per sample). Example of the data, which some packages store in the same object:
Cov:
      A      B      C      D      E
     19     10     18     14     20
     13      8     20     12     12
     12      6      9     16      7
     10      7      9     16      8
      8      5      6     12      8

M-values:
      A      B      C      D      E
   9.88   0.00   0.00   0.00   0.00
   7.93   0.00  10.00   3.00   0.00
   9.00   1.98   5.94  13.92   4.97
   9.00   2.94   3.96   8.96   4.96
   2.96   4.00   1.98   4.92   4.96

As you can see the coverage of each position is not the same for each sample.
However, I haven't found any specific method for calculating a PCA taking into account the coverage: there doesn't seem to be a method that normalizes the M-values by the coverage of each position.
I checked the distribution of coverage and 91% is covered by less than 14 reads, but some positions are covered by more than 300 reads.
What is the best way to normalize the position by coverage and plot the PCA ?

This is because of the statistical issue of Eigen value calculations - transformation are generally not welcome. In the changing age of unsupervised learning the concept of PCA is changing admittedly. You can superimpose coverage post PCA to assess with it is affecting the output, i.e. if it starts forming clusters there's a problem and your approach is justified.

I have experience in the analysis of EM-seq sequencing datasets but I honestly never really used M-values instead of classic fractions of methylation. One common approach is too remove highly covered positions (outliers) as it is often the results of sequencing or mapping errors. If you have a average coverage of ~ 14X, I think a threshold at a 100X is good. On the detection of methylation variation accross sample conditions, you won't need to deal with coverage variations a thes DSS soft (for DMC and DMR detection) works with raw counts. Not sure if it helps!

@PaulEndymion I also used the beta-values for a PCA (it revealed two other sample outliers). Interesting point about the mapping and the errors, I haven't actually mapped it so I don't know how well it was done. Are you suggesting that I don't need to perform PCA because they are already taken into account? I will look if removing those highly covered positions changes anything. Thanks!

Yes mapping errors on high repetition/GC% regions is a common bias, especially with short reads. These regions have the potential to attract reads with sequencing errors. We would do a PCA after the preprocessing of raw counts from "stranded_CpG_reports" files. Preprocessing included normalization of counts and removing of known SNP falling on "G" or "C" positions in "CG" motifs as well as higly covered positions. However as @M__ said if you found clusters of samples with PCA you can't explain with your samples experiment conditions, you might have a problem with these samples.

Just to be clear on the PCA. You can perform PCA but more as part of an exploratory/descriptive analysis, I don't remember reading about differential methylation studies using only PCA for DMR or DMC detection (even though DNA methylation studies are tedious so there is room for new approaches).


text:
Is it possible to create SMILES representation out of DNA/RNA-Oligonucleotides (~50 nts)?
I want to use physico chemical descriptors of single stranded DNA/RNA-Oligonucleotides in order to model certain properties. To generate these descriptors i need a computer understandable representation of these oligonucleotides, e.g. the SMILES representation. So far i could not find of any method to generate these SMILES representations out of the primary structure of these oligonucleotides.
I am wondering if someone knows of a software, that takes the primary sequence of a single stranded DNA-Oligonucleotide with e.g. 50 nts as input and gives the SMILES representation as output.
So far i could find nothing that would accomplish this task, although, at least to me, it seems like a straight forward task.
I will appreciate every constructive thought on that.
Thanks!
EDIT: Add some form of input and output file
As requested here an example of input and output pair:
input should be primary sequence:
TCA
output should be the corresponding SMILES string:
O=C(C(C)=CN1[C@@]2([H])C([H])([H])C@@([H])C@([H])O2)NC1=O
This example was generated with chemdraw 20.1.1.
For my use case i need a software, that convert many (like 50k - 1000k) primary sequences of single stranded DNA/RNA with a length of about 50 nts into the corresponding SMILES string. Ideally would be some sort of python package but anything else would help also a lot!

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.

Can you please provide an example input file and example output file?

It looks like various chemical software toolkits can export SMILES, for example OpenBabel, using the babel command line tool. I am not readily seeing the best way to read e.g. a FASTA, though. Another option would be ChemDraw.
This researchgate thread has a lot of information on peptide->SMILES conversion, I don't see nucleotides though.
One option, though it would be a little finicky, would be to look up the nucleotides SMILES strings on PubChem and then write a script to string them together (with a monophosphate) based on an input sequence. It seems like there is some related tooling already in rdkit, but a convenient API doesn't seem to exist yet. The developer does have one suggestion for peptides, which would translate into something like this for nucleotides:
from rdkit import Chem

nt_mol = Chem.MolFromFASTA("ACATATAT", flavor=9)  # flavor param notes as DNA
Chem.MolToSMILES(nt_mol)

To do this across many nucleotides you could hook it up to biopython to iterate across a FASTA file.
The nucleotide flavor option exists, but there is a note that only peptide is supported. I don't know whether that means the feature is experimental or is not implemented.


text:
Is there a software that extract molecular descriptors of single stranded RNAs secondary structure from dot bracket notation in nucleotide resolution?
I want to generate descriptors from dot bracket notations of single stranded RNA/DNA secondary structures. So far i have come across this paper:https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0097696#s5
In that, there is given a list of descriptors derived from the secondary structure of single stranded RNA but it is not explained how they derived these descriptors. Is there a Software out there, that would accomplish that task?
I wrote a script myself, to extract loop counts etc. from dot bracket notation but it would rather tedious to write a script all by myself for all descriptors possible.
I would appreciate every constructive thought on that.
Thanks!
EDIT Summary of paper for relevance
In this paper the secondary structures of single stranded RNA is parsed into many descriptors, whereby the secondary structure is given as dot bracket notation and where obtained by mfold. THis paper is nothing more than an example of the variety of descriptors that are describing the secondary structure of single stranded RNA/DNA present in dot bracket notation.
My question now is: Is there a general accepted software, that takes as input the dot bracket notation of an single stranded RND/DNA Oligonucleotide and returns many different descriptors? Such alike rdkit would do it for small organic molecules when given the SMILES string as input.
I hope this clarified my need.

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

Essentially, yes. Could you kindly summarise the PLoS article, for the information of relevance. Reading off-site articles to understand a question is not encouraged.

The authors use

Kapetanovic IM (2008) Computer-aided drug discovery and development (CADDD): in silico-chemico-biological approach. Chem Biol Interact 171: 165–176.

This generates the descriptors, which are normally used in small molecule interactions and they obtained a moderate (in reality low for ML) correlation. That is the innovative part of the study. "Decision trees" used alone are not recommended because it will over tighten onto the data. The reason they used them was to maximise the fit - the risk is over tightening. When its says "the algorithm successfully ...", well thats using a Decision Tree, any other ML algorithm will be lower.
The very nice part of the paper is the combination of wet-lab and analytical work, that balance is often not logistically trivial to achieve.
What is not clear is they talk about 6 descriptors out of 60 and also 27 descriptors out of 60. Which is the one they use? I would suggest the paper was not rigorously reviewed/edited to permit this sort of confusion for the reader. It could be (suspect) they performed two different analyses, course grain (6) and fine grain (27), but I don't think that is described very well.
I will not critique further.
Summary The paper is commendable for using a small molecular physiochemical/ thermodynamic description algorithm onto RNA data. Thats cool. Wrapping it around a Decision Tree alone, I wouldn't agree with - nothing wrong with DT's but they can't be used in isolation.

QSAR Model

In order to find the mathematical relationship between aptamers’ structure and activity in the training set, two functions of MASS packages of R program were used (R 2.14.1 2011 The R Foundation). Firstly, a stepwise-selection with the function “stepAIC” (direction = both) focused on the most relevant 27 descriptors out of the 60 that can contribute to the prediction of the relative activity. Then, a second step with multiple-linear-regression function “lm” was applied on the selected descriptors to generate appropriate coefficients for the linear correlation equation. The quality of the output was evaluated according to the p-value of the null hypothesis, and multiple R2 indicated the correlation between actual and predicted values.

QSAR

Sixty descriptors quantifying structural and thermodynamic properties of the aptamers were produced and gradually the most relevant descriptors were selected to form a six descriptor- based algorithm, which can predict the binding of aptamers to influenza virus and thus anticipate its anti-viral activity. The algorithm successfully predicted the binding affinity in the training set of aptamers (R2 = 0.702) and the was validated with a test set (R2 = 0.66); moreover, the model’s sensitivity and specificity in selecting aptamers with enhanced relative binding were high, 89% and 87%, respectively.

Personally, I don't think that is the question the OP is asking - I think they are asking about a different technique in RNA secondary structure analysis to the technique authors performed here. That however is simply a guess. It's an interesting paper, there is no question about that - it's not how I personally would do it - but  it's interesting. The issue is the thermodynamics has been rigorously worked out historically for RNA folding - it's a large body of work - and I'm sure the one or two reviewers would have pointed this out. It's complicated and if the authors are using their predictions to recycle back to the wet-lab .... thats what it's all about.

Normally, folding is calculated using graph theory ... FYI.

From there comments .... whats the program ... its Kapetanovic (2008) as cited above, which will be a classic small molecular binding algorithm the authors called QSAR.
In cases there's any misunderstanding - we wouldn't use small molecule interaction based calculations (MD) to understand RNA folding/binding because there's an alternative body of theory for this specific calculation, partly thermodynamics, partly a separate factor. ML is heavily used but not via an MD style calculation - which is what the authors are doing. Thus if you use the method there is a risk of criticism.
It's complicated - I understand the authors perspective, they are saying classical style hybridisation its just not that simple. It might be complicated but hybridisation shouldn't be ignored as the empirical assumption.

♦ Thanks for the comment! Do you know of a software that generates such descriptors when given the dot bracket notation? Ideally would be some sort of python package. Many Descriptors in nucleotide resolution are the product of more or less simple string parsing of the dot bracket notation and i would rather not  write these parsings all by myself.

Again @Gojih the authors cite Kapetanovic IM (2008) Computer-aided drug discovery and development (CADDD): in silico-chemico-biological approach. Chem Biol Interact 171: 165–176. This is the basis of their QSAR analysis. Its basically MD for nucleic acids

@Gojih answered in comments.


text:
Error while using SummarizedExperiment() in R
I'm tryig to perform RKM normalization on raw counts for RNA-Seq Data:
counts_data <- exprs181517[,-1]
#Setting title col to rownames for pheno data
rownames(pheno.data181517) <- NULL
rownames(counts_data) <- exprs181517[,1]
phenoColData <- tibble::column_to_rownames(pheno.data181517, "title")
colnames(counts_data) <- rownames(phenoColData)
#Defining the design column
designCol <- as.data.frame(phenoColData$characteristics_ch1)
rownames(designCol) <- rownames(phenoColData)

#Creating a summarized experiment
se <- SummarizedExperiment(list(counts=counts_data), colData=phenoColData)

DESeqDataSetFromMatrix(se, design ~ designCol)

but I keep getting this error:
Error in h(simpleError(msg, call)) :
  error in evaluating the argument 'x' in selecting a method for function 'as.matrix':
    no method for coercing this S4 class to a vector

DESeqDataSetFromMatrix but a SummarizedExperiment is not a matrix. Use the function DESeqDataSet() instead and make sure the design argument is followed by a =.


text:
How to Recalculate Allele Frequency from a VCF File?
I have a VCF file with 200 samples (mitochondrial genome of Plasmodium falciparum). Here is a pic to take a look at:

And a few relevant lines from the actual file:
##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed">
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed">
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  FP0008-C   FP0009-C
Pf_M76611       19      .       A       T       33196.8 MissingVQSLOD;Mitochondrion     AC=9;AF=0.0003019;AN=29816;ANN=T|intragenic_variant|MODIFIER|PF3D7_MIT00100|PF3D7_MIT00100|gene_variant|PF3D7_MIT00100|||n.19A>T||||||;AS_InbreedingCoeff=0.0437;AS_QD=20.96;BaseQRankSum=0.141;CDS;DP=559138;ExcessHet=0.0026;FS=2.892;InbreedingCoeff=0.52;MLEAC=8;MLEAF=0.0002683;MQ=60;MQRankSum=0;QD=20.96;ReadPosRankSum=0.856;RegionType=Mitochondrion;SOR=0.789;set=genotypegvcfs  GT:AD:DP:GQ:PL  0/0:52,0:52:99:0,120,1800  0/0:54,0:54:99:0,120,1800

A small subset of samples I am working on can be found here:
https://www.mediafire.com/file/6dkct6guk5zx83m/sample.vcf/file
I am trying to calculate allele frequency (AF) for each variant in the dataset. The population is assumed to be all the samples in the file.
In other words, even though the VCF file contains AF statistics from population data, I would like to know how to recalculate AF only from the VCF file. I believe I need to know how AC and AN are calculated first.
AF is calculated this way: allele count (AC) / allele number (AN)
I tried to calculate them from the information from GT (counting 0/0, 1/1 and 0/1 in all the 200 samples), but then I learned that this calculation is not suitable for mithocondrial data.
All lines in the VCF file contain an AF statistic:
$ grep -v 'AF=' sample.vcf | grep -v '^#'
<No output>

The information is provided under INFO column, but how were they calculated? If I do not have AC and AN to begin with, what should I do?
Is there a way to calculate AC and AN from the other information provided in the dataset? such as the data mentioned in the FORMAT column?
GT: genotype
AD: unfiltered allele depth
DP : read depth at this position for this sample (Integer)
GQ : conditional genotype quality, encoded as a phred quality
PL : the phred-scaled genotype likelihoods rounded to the closest integer
Thanks

Can you give us some more context? What is AC? Allele count in the population? From GnomAD, maybe? What tool produced this file?

@terdon Thanks for the reply. I believe that is what AC means. I have not seen any other definition.

It sounds like your main confusion is how to calculate AF. I'll answer that right away, but I'll also suggest a Python library that you can use to do this and much more.
You're correct, the AF is defined as the number of non-missing alternate alleles (AC) divided the number of non-missing alleles (AN). Calculating these are straightforward. Map each GT to its alternate alleles:
0/0 => 0
0/1 => 1
1/1 => 2
./0 => 0
./1 => 1
./. => 0

Now, compute the sum. That is the AC. For the AN, map each GT to its number of non-missing alleles:
0/0 => 2
0/1 => 2
1/1 => 2
./0 => 1
./1 => 1
./. => 0

Now, compute the sum. That is the AN. Now calculate AF as
AC / AN

All of this is kind of annoying to do though! I suggest you use Hail (disclaimer: I am a Hail developer).
import hail as hl

# Hail is not yet aware of Plasmodium falciparum reference genome
mt = hl.import_vcf('sample.vcf', reference_genome=None)

# let's peek at our data
mt.show(n_rows=10, n_cols=5)

# Let's get some information about the GT specifically:
mt.GT.summarize()

# We can use Hail's built in variant_qc method to calculate AF for us:
mt = hl.variant_qc(mt)

# Let's take a look at that
mt.variant_qc.show()

# We can overwrite the old AF field and export back to a VCF file if necessary
mt = mt.annotate_rows(
    # AF[1] is the frequency of the first alternate allele
    # Your dataset appears to have some multi-allelics. The other
    # alternate allele frequencies are available in AF[2], ...
    info = mt.info.annotate(
        AF = mt.variant_qc.AF[1]
    )
)
hl.export_vcf(mt, 'sample-updated.vcf.bgz')

This comment doesn't make sense to me:

I learned that this calculation is not suitable for mithocondrial data.

The definition of AF, AC, and AN doesn't change depending on the type of config (autosome, sex chromosome, mitochondrial). Perhaps the issue is ploidy encoding? A GT in a VCF can be:

1-ploidy, e.g. 0 (reference), 1 (first alternate), 2 (second alternate)
2-ploidy, e.g. 0/0 (homozygous reference), 0/1 (heterozygous), 1/1 (homozygous first alternate)
etc.

Sometimes people incorrectly encode a 1-ploidy config (for example the Y chromosome) as: 0/0 and 1/1 instead of 0 and 1. That's not your fault though and it doesn't mean the definition of AF, AC, and AN change. It means someone gave you malformed data.

Thanks for the reply. You totally understood my question. Reagrding the part you said that does not make sense, what i meant was that cells can have multiple mts. From the info i gathered so far, due to this reason the simple equation of AF = AC / AN would not be appropriate for mt genome meaning using GT alone will not return the right result. The calculations used in softwares such as vcftools are only for non mt genome. I was trying to do it in Python anyway. I am glad you introduced this library. Does it work on Windows too?

We don't regularly test against windows, but, in principle, it should work. Let me know if you run into bugs and we'll fix them!


text:
How is the "canonical" version (`_entity_poly.pdbx_seq_one_letter_code`) obtained in the PDB?
I encountered the dichotomy in the context of PDBx/mmcif files, say, 6OSQ: each chain has a
_entity_poly.pdbx_seq_one_letter_code but also _entity_poly.pdbx_seq_one_letter_code_can which i assume stands for "canonical". The raw and the canonical sequences.
The difference seems to be that all the ligands, hetatoms and other junk ((PSU), (MP4), (MG4)) have been removed from the the raw sequence to make it canonical i.e. it can now be packaged as fasta.
My question is where does one find documentation for this process? Is making it amenable to fasta tools the primary motivation? What kind of hetero atoms can occur there? Where is the logic written down for replacing them or deleting them when deriving the canonical sequence? Is the job of the structure's depositor or the PDB? Is there a tool for this?
Below is chain 1 of PDB6OSQ. First is the raw seq, then is the canonical one.
1 polyribonucleotide no yes
;GGUUAAGCGACUAAGCGUACACGGUGGAUGCCCUGGCAGUCAGAGGCGAUGAAGGACGUGCUAAUCUGCGAUAAGCGUCG
GUAAGGUGAUAUGAACCGUUAUAACCGGCGAUUUCCGAAUGGGGAAACCCAGUGUGUUUCGACACACUAUCAUUAACUGA
AUCCAUAGGUUAAUGAGGCGAACCGGGGGAACUGAAACAUCUAAGUACCCCGAGGAAAAGAAAUCAACCGAGAUUCCCCC
AGUAGCGGCGAGCGAACGGGGAGCAGCCCAGAGCCUGAAUCAGUGUGUGUGUUAGUGGAAGCGUCUGGAAAGGCGCGCGA
UACAGGGUGACAGCCCCGUACACAAAAAUGCACAUGCUGUGAGCUCGAUGAGUAGGGCGGGACACGUGGUAUCCUGUCUG
AAUAUGGGGGGACCAUCCUCCAAGGCUAAAUACUCCUGACUGACCGAUAGUGAACCAGUACCGUGAGGGAAAGGCGAAAA
GAACCCCGGCGAGGGGAGUGAAAAAGAACCUGAAACCGUGUACGUACAAGCAGUGGGAGCACGCUUAGGCGUGUGACUGC
GUACCUUUUGUAUAAUGGGUCAGCGACUUAUAUUCUGUAGCAAGGUUAACCGAAUAGGGGAGCCGAAGGGAAACCGAGUC
UUAACUGGGCGUUAAGUUGCAGGGUAUAGACCCGAAACCCGGUGAUCUAGCCAUGGGCAGGUUGAAGGUUGGGUAACACU
AACUGGAGGACCGAACCGACUAAU(1MG)(PSU)(5MU)GAAAAAUUAGCGGAUGACUUGUGGCUGGGGGUGAAAGGCCA
AUCAAACCGGGAGAUAGCUGGUUCUCCCCGAAAGCUAUUUAGGUAGCGCCUCGUGAAUUCAUCUCCGGGGGUAGAGCACU
GUUUCGGCAAGGGGGUCAACCCGACUUACCAACCCGAUGCAAACUGCGAAUACCGGAGAAUGUUAUCACGGGAGACACAC
GGCGGG(PSU)GCUAACGUCCGUCGUGAAGAGGGAAACAACCCAGACCGCCAGCUAAGGUCCCAAAGUCAUGGUUAAGUG
GGAAACGAUGUGGGAAGGCCCAGACAGCCAGGAUGUUGGCUUAGAAGCAGCCAUCAUUUAAAGAAAGCGUAAUAGCUCAC
UGGUCGAGUCGGCCUGCGCGGAAGAUGUAACGGGGCUAAACCAUGCACCGAAGCUGCGGCAGCGACGCUUAUGCGUUGUU
GGGUAGGGGAGCGUUCUGUAAGCCUGCGAAGGUGUGCUGUGAGGCAUGCUGGAGGUAUCAGAAGUGCGAAUGCUGACAUA
AGUAACGAUAAAGCGGGUGAAAAGCCCGCUCGCCGGAAGACCAAGGGUUCCUGUCCAACGUUAAUCGGGGCAGGGUGAGU
CGACCCCUAAGGCGAGGCCGAAAGGCGUAGUCGAUGGGAAACAGGUUAAUAUUCCUGUACUUGGUGUUACUGCGAAGGGG
GGACGGAGAAGGCUAUGUUGGCCGGGCGACGGUUGUCCCGGUUUAAGCGUGUAGGCUGGUUUUCCAGGCAAAUCCGGAAA
AUCAAGGCUGAGGCGUGAUGACGAGGCACUACGGUGCUGAAGCAACAAAUGCCCUGCUUCCAGGAAAAGCCUCUAAGCAU
CAGGUAACAUCAAAUCGUACCCCAAACCGACAC(6MZ)GGUGGUCAGGUAGAGAAUACCAAGGCGCUUGAGAGAACUCGG
GUGAAGGAACUAGGCAAAAUGGUGCCGUAACUUCGGGAGAAGGCACGCUGAUAUGUAGGUGAGGUCCCUCGCGGAUGGAG
CUGAAAUCAGUCGAAGAUACCAGCUGGCUGCAACUGUUUAUUAAAAACACAGCACUGUGCAAACACGAAAGUGGACGUAU
ACGGUGUGACGCCU(2MG)CCCGGUGCCGGAAGGUUAAUUGAUGGGGUUAGCGCAAGCGAAGCUCUUGAUCGAAGCCCCG
GUAAACGGCGGCCG(PSU)AAC(3TD)A(PSU)AACGGUCCUAAGGUAGCGAAA(5MU)UCCUUGUCGGGUAAGUUCCGA
C(5MC)UGCACGAAUGGCGUAAUGAUGGCCAGGCUGUCUCCACCCGAGACUCAGUGAAAUUGAACUCGCUGUG(6MZ)AG
AUGCAGUGUACCCGCGGCAAGACGGAAAGACCCCGU(G7M)AACCUUUACUAUAGCUUGACACUGAACAUUGAGCCUUGA
UGUGUAGGAUAGGUGGGAGGCUUUGAAGUGUGGACGCCAGUCUGCAUGGAGCCGACCUUGAAAUACCACCCUUUAAUGUU
UGAUGUUCUAACGUUGACCCGUAAUCCGGGUUGCGGACAGUGUCUGGUGGGUAGUUUGACUG(OMG)GGCGGUCUCCUCC
UAAAGAGUAACGGAGGAGCACGAAGGUUGGCUAAUCCUGGUCGGACAUCAGGAGGUUAGUGCAAUGGCAUAAGCCAGCUU
GACUGCGAGCGUGACGGCGCGAGCAGGUGCGAAAGCAGGUCAUAGUGAUCCGGUGGUUCUGAAUGGAAGGGCCAUCGCUC
AACGGAUAAAAGGUACUCCG(2MG)GGAUAACAGGC(PSU)GAUACCGCCCAAGAGUUCAUAUCGACGGCGGUGUUUGGC
A(OMC)CUCG(2MA)(PSU)GUCGGCUCAUCACAUCCUGGGGCUGAAGUAGGUCCCAAGGGUAUGGC(OMU)GUUCGCCA
UUUAAAGUGGUACGCGAGC(PSU)GGGUUUAGAACGUCGUGAGACAGU(PSU)CGGUCCCUAUCUGCCGUGGGCGCUGGA
GAACUGAGGGGGGCUGCUCCUAGUACGAGAGGACCGGAGUGGACGCAUCACUGGUGUUCGGGUUGUCAUGCCAAUGGCAC
UGCCCGGUAGCUAAAUGCGGAAGAGAUAAGUGCUGAAAGCAUCUAAGCACGAAACUUGCCCCGAGAUGAGUUCUCCCUGA
CCCUUUAAGGGUCCUGAAGGAACGUUGAAGACGACGACGUUGAUAGGCCGGGUGUGUAAGCGCAGCGAUGCGUUGAGCUA
ACCGGUACUAAUGAACCGUGAGGCUUAACCU
;
;GGUUAAGCGACUAAGCGUACACGGUGGAUGCCCUGGCAGUCAGAGGCGAUGAAGGACGUGCUAAUCUGCGAUAAGCGUCG
GUAAGGUGAUAUGAACCGUUAUAACCGGCGAUUUCCGAAUGGGGAAACCCAGUGUGUUUCGACACACUAUCAUUAACUGA
AUCCAUAGGUUAAUGAGGCGAACCGGGGGAACUGAAACAUCUAAGUACCCCGAGGAAAAGAAAUCAACCGAGAUUCCCCC
AGUAGCGGCGAGCGAACGGGGAGCAGCCCAGAGCCUGAAUCAGUGUGUGUGUUAGUGGAAGCGUCUGGAAAGGCGCGCGA
UACAGGGUGACAGCCCCGUACACAAAAAUGCACAUGCUGUGAGCUCGAUGAGUAGGGCGGGACACGUGGUAUCCUGUCUG
AAUAUGGGGGGACCAUCCUCCAAGGCUAAAUACUCCUGACUGACCGAUAGUGAACCAGUACCGUGAGGGAAAGGCGAAAA
GAACCCCGGCGAGGGGAGUGAAAAAGAACCUGAAACCGUGUACGUACAAGCAGUGGGAGCACGCUUAGGCGUGUGACUGC
GUACCUUUUGUAUAAUGGGUCAGCGACUUAUAUUCUGUAGCAAGGUUAACCGAAUAGGGGAGCCGAAGGGAAACCGAGUC
UUAACUGGGCGUUAAGUUGCAGGGUAUAGACCCGAAACCCGGUGAUCUAGCCAUGGGCAGGUUGAAGGUUGGGUAACACU
AACUGGAGGACCGAACCGACUAAUGUUGAAAAAUUAGCGGAUGACUUGUGGCUGGGGGUGAAAGGCCAAUCAAACCGGGA
GAUAGCUGGUUCUCCCCGAAAGCUAUUUAGGUAGCGCCUCGUGAAUUCAUCUCCGGGGGUAGAGCACUGUUUCGGCAAGG
GGGUCAACCCGACUUACCAACCCGAUGCAAACUGCGAAUACCGGAGAAUGUUAUCACGGGAGACACACGGCGGGUGCUAA
CGUCCGUCGUGAAGAGGGAAACAACCCAGACCGCCAGCUAAGGUCCCAAAGUCAUGGUUAAGUGGGAAACGAUGUGGGAA
GGCCCAGACAGCCAGGAUGUUGGCUUAGAAGCAGCCAUCAUUUAAAGAAAGCGUAAUAGCUCACUGGUCGAGUCGGCCUG
CGCGGAAGAUGUAACGGGGCUAAACCAUGCACCGAAGCUGCGGCAGCGACGCUUAUGCGUUGUUGGGUAGGGGAGCGUUC
UGUAAGCCUGCGAAGGUGUGCUGUGAGGCAUGCUGGAGGUAUCAGAAGUGCGAAUGCUGACAUAAGUAACGAUAAAGCGG
GUGAAAAGCCCGCUCGCCGGAAGACCAAGGGUUCCUGUCCAACGUUAAUCGGGGCAGGGUGAGUCGACCCCUAAGGCGAG
GCCGAAAGGCGUAGUCGAUGGGAAACAGGUUAAUAUUCCUGUACUUGGUGUUACUGCGAAGGGGGGACGGAGAAGGCUAU
GUUGGCCGGGCGACGGUUGUCCCGGUUUAAGCGUGUAGGCUGGUUUUCCAGGCAAAUCCGGAAAAUCAAGGCUGAGGCGU
GAUGACGAGGCACUACGGUGCUGAAGCAACAAAUGCCCUGCUUCCAGGAAAAGCCUCUAAGCAUCAGGUAACAUCAAAUC
GUACCCCAAACCGACACAGGUGGUCAGGUAGAGAAUACCAAGGCGCUUGAGAGAACUCGGGUGAAGGAACUAGGCAAAAU
GGUGCCGUAACUUCGGGAGAAGGCACGCUGAUAUGUAGGUGAGGUCCCUCGCGGAUGGAGCUGAAAUCAGUCGAAGAUAC
CAGCUGGCUGCAACUGUUUAUUAAAAACACAGCACUGUGCAAACACGAAAGUGGACGUAUACGGUGUGACGCCUGCCCGG
UGCCGGAAGGUUAAUUGAUGGGGUUAGCGCAAGCGAAGCUCUUGAUCGAAGCCCCGGUAAACGGCGGCCGUAACXAUAAC
GGUCCUAAGGUAGCGAAAUUCCUUGUCGGGUAAGUUCCGACCUGCACGAAUGGCGUAAUGAUGGCCAGGCUGUCUCCACC
CGAGACUCAGUGAAAUUGAACUCGCUGUGAAGAUGCAGUGUACCCGCGGCAAGACGGAAAGACCCCGUGAACCUUUACUA
UAGCUUGACACUGAACAUUGAGCCUUGAUGUGUAGGAUAGGUGGGAGGCUUUGAAGUGUGGACGCCAGUCUGCAUGGAGC
CGACCUUGAAAUACCACCCUUUAAUGUUUGAUGUUCUAACGUUGACCCGUAAUCCGGGUUGCGGACAGUGUCUGGUGGGU
AGUUUGACUGGGGCGGUCUCCUCCUAAAGAGUAACGGAGGAGCACGAAGGUUGGCUAAUCCUGGUCGGACAUCAGGAGGU
UAGUGCAAUGGCAUAAGCCAGCUUGACUGCGAGCGUGACGGCGCGAGCAGGUGCGAAAGCAGGUCAUAGUGAUCCGGUGG
UUCUGAAUGGAAGGGCCAUCGCUCAACGGAUAAAAGGUACUCCGGGGAUAACAGGCUGAUACCGCCCAAGAGUUCAUAUC
GACGGCGGUGUUUGGCACCUCGAUGUCGGCUCAUCACAUCCUGGGGCUGAAGUAGGUCCCAAGGGUAUGGCUGUUCGCCA
UUUAAAGUGGUACGCGAGCUGGGUUUAGAACGUCGUGAGACAGUUCGGUCCCUAUCUGCCGUGGGCGCUGGAGAACUGAG
GGGGGCUGCUCCUAGUACGAGAGGACCGGAGUGGACGCAUCACUGGUGUUCGGGUUGUCAUGCCAAUGGCACUGCCCGGU
AGCUAAAUGCGGAAGAGAUAAGUGCUGAAAGCAUCUAAGCACGAAACUUGCCCCGAGAUGAGUUCUCCCUGACCCUUUAA
GGGUCCUGAAGGAACGUUGAAGACGACGACGUUGAUAGGCCGGGUGUGUAAGCGCAGCGAUGCGUUGAGCUAACCGGUAC
UAAUGAACCGUGAGGCUUAACCU

You can check specification of pdbx_seq_one_letter_code and pdbx_seq_one_letter_code_can. These items are generated from the depositor-provided sequence on the PDB side. I suppose that the software used is MAXIT.
Note: these are not PDB files. They are distributed by the PDB organization, but the file format is mmCIF (or, more formally, PDBx/mmCIF).

Thanks Marcin, love the gemmi package. You are right, i was talking about mmcif imprecisely.
I did find those references, but they seem a little incomplete at least.. For example they don't describe this (PSU) element which i found in literally the first mmcif file on my computer that i opened. I know now that (PSU) is pseudo-uracil because i stumbled across some .txt doc file, but i can't find it anymore.

So, i was little myopic when reading the PDB description the first few times, but roughly what i was looking for is this:

the non-canonical sequence(pdbx_seq_one_letter_code) makes explicit the ligands present in the given chain.
the full dictionary of these is called the CCD and contains ~20000 entries, is curated and updated weekly by the PDB. That's basically their database (glorified text file over ftp) against which they check whether a given ligand has already been recorded.
the process of chaffing them out of the EM dough as well as schema and more context canbe found here: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4243272/pdf/bau116.pdf

... the non-canonical sequence(pdbx_seq_one_letter_code) makes explicit the ligands present in the given chain.... now you have to explain how ligands are assigned to relevant chain A instead of equally close one chain BB ?


text:
Is there a smaller representative version of refseq with fewer than 30K files in it?
Does any one have a condensed NCBI reference seq database with all the gbff files, which covers mostly all the bacteria and archaea. The current NCBI ReqSeq contains more than 30,000 gbff files. I am looking for a much more curated refseq with less than 30K files in it.

Please add more context to this question to explain why you need a reduced refseq dataset. There are lots of different reasons for this, and those reasons will change how the reduction is carried out.

proGenomes includes a collection of genomes that are representative of different types - https://progenomes.embl.de/download.cgi
Or Genome Updater could be used to filter on unique taxonomy IDs and genome completeness etc.
https://github.com/pirovc/genome_updater


text:
Is the potential skewing of FDR and ranking of gene set analysis with gene sets with hierarchical structure a concern?
When performing GSEA and many other gene set analysis, when gene sets with hierarchy such as KEGG collection or Reactome collection, the FDR could be skewed sometimes, and the ranking of gene sets by any metrics could be confusing because multiple subsets belonging to the same set would have similar ranking, pushing other gene sets down the list. For example, in reactome axon guidance (R-HSA-422475.6) contains semaphorin interaction (R-HSA-373755.1), which contains Sema4d in semaphorin signaling (R-HSA-400685.2). The enrichment of Sema4d gene set usually (though not necessarily) entails the enrichment of the other gene sets above it. That might lead to an over-adjustment of p-value if I am not mistaken.
I am just wondering if there are ways to mitigate it, or it is not a major concern that I should have?

This is actually a known problem in the correction of nominal p-values in over-representation analyses that use the Gene Ontology categories (commonly but misleadingly called "GO enrichment" in many publications).
Most currently used methods for p-value adjustment (like BH) assume that p-values are independent, which is not true for a case like the one you describe - a set contained by another set is basically by definition not independent. Even if you didn't have a hierarchical structure, you would still have many overlaps across all your genesets, which would still mean that you cannot assume independence. So in both cases, but I think especially in your case, a BH/FWER correction would be statistically wrong (and possibly overly conservative).
There are a few ways to avoid this: one way would be to prune the directed acyclic graph of your genesets/categories, retaining only the categories that are enriched and have no enriched children.
In your case, for instance, if "axon guidance" and  "semaphorin interaction" are significant but "Sema4d in semaphorin signaling" is not, you would stop at "semaphorin signaling" removing the result for "axon guidance".
This is kind of a strong correction, which is guaranteed to lower the multiple testing correction burden (unless you have no significant genesets).
There are however more refined ways to address this:
One of them is to remove from set A (parent) any genes that belong to a significantly over-represented/enriched set A' (child). This method is termed the elim method and was described in this publication.
The same publication discusses another method called weight which assigns a weight to a significant category based on the significance of its children.
Another method is implemented in the gProfiler tool and introduced in this 2007 publication, which is the g:SCS method. Namely it calculates p-value thresholds based on empirical, randomized p-values for sets of different sizes. According to the authors it is more conservative than BH but less than Bonferroni.
The only issue with these approaches is that there are, to my knowledge, no ways to apply them to your data independently, i.e. they are coded within the respective tools and do not exist as stand-alone packages for multiple test correction. I hope I'm wrong and someone developed these versions, but for the time being the only thing you can do is to try to implement your own versions of these methods.
Another relevant, more recent publication is the SetRank one, which describes a method to do GSEA with a p-value calculation that accounts for overlaps, and has its own package on CRAN. If you prefer to apply a correction to your own results (like the ones you get from fgsea), however, I do not know whether this package has methods you can use or recycle.

Great answer! Is there a recommended standard package in R to trim DAGs without that one codes that up yourself?

@ATpoint IMHO it depends on the type of test. For Over-Representation Analysis (ORA), which is the one that normally uses GO genesets, packages such as topGO (Bioconductor) and gProfiler2 (CRAN) already implement trimming. For GSEA - what the user was asking about - SetRank may be one choice, but I have not used it personally. I am partial to fgsea for speed, so I did code my own pruning functions some time ago. I believe they are a little crude as they remove entire categories based on the DAG, but don't really address overlaps.

Thank you! Both answers provide some very good insight in gene set analysis. Trimming based on DAGs seems very straightforward and easy to do given database such as Reactome provides them in various format already. While searching, I also found a package called GeneSetCluster which takes GSA output and correct the results for overlapping genes between gene sets. I haven't tried it yet, but this could be an addition the DAG trimming and SetRank.

Statisticians might kill me for it but I've always filtered long and redundant geneset lists to only contain what remotely makes biological sense in my setup. If I am investigating leukocyte transcriptomes I do not care about terms related to filtration processes in the kidney or signal transduction in the hypothalamus. I think the redundancy of common pathway databases in terms of genes contributing to many terms justifiers to eliminate certain terms and choose those that make sense and are not overly large (for example not more than 500 genes). Top-of-hierarchy terms like 'metabolism' eith more than 1000 genes are anyway unspecific and hard to interpret. That then reduces the multiple testing burden quite a lot.


text:
How to handle output from a python script that gives different outputs depending upon the condition - Nextflow
I have created three Nextflow processes: One checks the strandedness of the fastq file and outputs it to a .txt file. Second one executes a python script with different commands depending upon the condition. And the third one aligns them.
/*
 * the `CHECK_STRANDEDNESS` process takes SAM file
 * and a BED file
 * and outputs a a txt file containing information about strandedness.
 */

process CHECK_STRANDEDNESS {

    input:
    path "${pair_id}.sam"
    path "${gtf.baseName}.bed"

    output:
    path "${pair_id}_strand_check.txt"

    script:
    """
    infer_experiment.py -i "${pair_id}.sam" -r "${gtf.baseName}.bed" > "${pair_id}_strand_check.txt"
    """
}

/*
 * the `FINDTAIL` process takes "${pair_id}_strand_check.txt" file
 * and executes commands according to the reads strandedness.
 * The output generated is FASTA files.
 */

process FINDTAIL {

    input:
    path "${pair_id}_strand_check.txt"
    path tuple val(pair_id), path(reads)

    output:

    script:
    """
    #!/usr/bin/env python3

    import subprocess
    import pandas as pd
    import os

    result = pd.read_csv("${pair_id}_strand_check.txt", sep="\r\n", header=None, engine='python')

    failed = float(result.iloc[1,0].replace('Fraction of reads failed to determine: ', ''))
    fwd = float(result.iloc[2,0].replace('Fraction of reads explained by "1++,1--,2+-,2-+": ', ''))
    rev = float(result.iloc[3,0].replace('Fraction of reads explained by "1+-,1-+,2++,2--": ', ''))
    fwd_percent = fwd/(fwd+rev)
    rev_percent = rev/(fwd+rev)

    if float(result.iloc[1,0].replace('Fraction of reads failed to determine: ', '')) > 0.50:
        cmd1 = "findtail_v1.01 --input_file ${reads[0]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_1.fasta
        cmd2 = "findtail_v1.01 --input_file ${reads[1]}  --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_2.fasta
        print(cmd1)
        subprocess.call(cmd1, shell=True)
        print(cmd2)
        subprocess.call(cmd2, shell=True)

        if fwd_percent > 0.9:
            #Over 90% of reads explained by "1++,1--,2+-,2-+
            #Data is likely FR/fr-secondstrand
            cmd1 = "findtail_v1.01 --input_file ${reads[0]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type pr --d > tail_${pair_id}_1.fasta
            cmd2 = "findtail_v1.01 --input_file ${reads[1]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type pr --d > tail_${pair_id}_2.fasta
            print(cmd1)
            subprocess.call(cmd1, shell=True)
            print(cmd2)
            subprocess.call(cmd2, shell=True)

        elif rev_percent > 0.9:
            # Over 90% of reads explained by "1+-,1-+,2++,2--
            # Data is likely RF/fr-firststrand
            cmd1 = "findtail_v1.01 --input_file ${reads[0]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_1.fasta
            cmd2 = "findtail_v1.01 --input_file ${reads[1]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_2.fasta
            print(cmd1)
            subprocess.call(cmd1, shell=True)
            print(cmd2)
            subprocess.call(cmd2, shell=True)

        elif max(fwd_percent, rev_percent) < 0.6:
            #Under 60% of reads explained by one direction
            #Data is likely unstranded
            cmd1 = "findtail_v1.01 --input_file ${reads[0]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_1.sl.fasta
            cmd2 = "findtail_v1.01 --input_file ${reads[0]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type pr --d > tail_${pair_id}_1.pr.fasta
            cmd3 = "findtail_v1.01 --input_file ${reads[1]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_2.sl.fasta
            cmd4 = "findtail_v1.01 --input_file ${reads[1]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type pr --d > tail_${pair_id}_2.pr.fasta
            print(cmd1)
            subprocess.call(cmd1, shell=True)
            print(cmd2)
            subprocess.call(cmd2, shell=True)
            print(cmd3)
            subprocess.call(cmd3, shell=True)
            print(cmd4)
            subprocess.call(cmd4, shell=True)

        else:
            #if strand couldnt be detected
            # we assume it is first-stranded
            # since most illumina reads are
            cmd1 = "findtail_v1.01 --input_file ${reads[0]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_1.fasta
            cmd2 = "findtail_v1.01 --input_file ${reads[1]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_2.fasta
            print(cmd1)
            subprocess.call(cmd1, shell=True)
            print(cmd2)
            subprocess.call(cmd2, shell=True)
        """
}

/*
 * define the `ALIGN` process that outputs unmapped files in FASTA format
 * It will take the FASTA files creates by the `FINDTAIL` process
 * and try to align them INDIVIDUALLY.
 */

process ALIGN {

    input:
    path index
    path "tail_${pair_id}*.fasta"

    output:
    path "unmapped_tail_${pair_id}*.fasta"

    script:
    """
    bowtie2 -x $index -f "tail_${pair_id}*.fasta" --un "unmapped_tail_${pair_id}*.fasta"
    """
}

After creating a "${pair_id}_strand_check.txt" file, I input it to the FINDTAIL process that finds reads with polyadenylated ends. It takes fastq files individually, so I have to run it individually for each file.
If fastq file pair is first stranded or second stranded, I executed findtail twice,once for R1 and once for R2, which creates two fasta files. But if it is unstranded, I have to run findtail four times, twice for each file. Thus, depending on the paired-end reads' strandedness, I have different number of inputs. How do I define the output of this process?
In the next process I align these files individually to get unmapped reads.
This issue has been eating my brain away for 5 days. So please let me know if something is unclear.
Best wishes!

From the (now removed) edits, the updated path glob pattern could just be: path("tail_${pair_id}_{1,2}.{pr,sl}.fasta"). If it's still not working, check to see what files are actually being created. They might be being created in another directory or not at all.

I ran it the same thing again and it worked somehow. So  gathered it was fault on my end and remove the edit.

You should be able to just use a glob pattern that accounts for all possible output files, regardless of the specified conditions. As long as a single file matches the glob pattern when the process completes (successfully), Nextflow won't complain. I think this should suffice, but please let me know in the comments if it doesn't:
output:
tuple val(pair_id), path("tail_${pair_id}_{1,2}{,.pr,.sl}.fasta")

Thank you so much! Merry Christmas and Happy Holidays!


text:
How to manually trim FASTA file sequences with the information provided in the header and store it into a new FASTA file? - Python
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!

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?

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

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?

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\w+),' r'start=(?P\d+),' r'end=(?P\d+),' r'length=(?P\d+),' r'identity=(?P\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?

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()

Its worth point out that re is notoriously slow and over large files parallelisation is needed across the Biopython object.

Hello Steve. Thank you so much for you reply but I tried running your code and it gave the following error: AttributeError: 'NoneType' object has no attribute 'groupdict'

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.

Hello @Steve, I have edited the question. If you could please take a look, thank you.

@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.

@pubsurfted There is indeed, using the "fasta-2line" format.

Thank you so so much Steve! Also, the "fasta-2line" format is raising the error: ValueError: Unknown format 'fasta-2line'

It ran. I had to update the biopython package from 1.67 to 1.80.