text:
Remove low quality reads
I want to remove reads from FASTQ file that contain homopolymers > 10bp and remove reads with <35 average quality score across the entire read.
To remove homopolymers > 10bp, I tried this on a Linux machine, but it only removes the sequence line:
zcat file.fastq.gz |
  awk '!/A{10,}/&&!/C{10,}/&&!/G{10,}/&&!/T{10,}/ {print}' cleaned_file.fastq

Please [edit] your question and explain what you need in more detail. What do you mean by "less than the 35 quality score"? Do you want to remove all reads that have even a single base with <35 score? Do you want to calculate the average quality across the entire read and remove reads with <35 average quality score? Something else? Why do you want to remove these from the fastq file instead of handling it in the aligned bam or as part of your variant calling process? How are you defining homopolymers here? Do you just want to remove reads if they have >10 consecutive identical bases?

Use an off the shelf tool for read preprocessing.
Here is one:
./fastq_qual_trimmer -i test.fq -m 35 -H 10

That one does exactly what you want but seems a little old/unmaintained, so here is another that does more or less the same thing:
fastq-mcf --qual-mean 35 --homopolymer-pct {X} adapters.fa reads.fq

where {X} is 10 / read length, adapters.fa is an adapter file (which I believe can be empty or filled with dummy sequences).
You could also use a library like biopython or dnaio to write a quick script to do this, but it hardly seems worth it.

Thanks for your help. Sorry I didn't write in detail but i need to calculate the average quality over the whole read and remove the reads with average quality score <35. As I understand it, these tools can't handle it.

@n_n sorry, used wrong options. Both of these tools can handle it across the read. Updated answer.

Thank you so much for taking the time to help me. Your assistance is greatly appreciated.

@user_01 if this answer has solved your problem could you very kindly please mark this as accepted? It would be a small way of showing gratitude and helps the site in several ways.

Just one more question I tried fastq-mcf as you suggest for my fastq files but i got this Error: number of input files must match number of '-o' output files. May you help me to understand whats wrong with code.
   read_length=$(head -n 2 ${file} | tail -n 1 | wc -m)
   homopolymer_pct=$((8/${read_length}))
   fastq-mcf ${file} adapters.fa --qual-mean 35 --homopolymer-pct  ${homopolymer_pct} -o filtered/${file}
done```

@user_01 I don't have the time right now to extensively troubleshoot but i'd guess that something is wrong with the way that you're using {$file}, or with the files you are passing, or with the existence of directory filtered/. If you print out {$file} or some of the other values you might see something.


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"

source: Stack Exchange

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

sort: votes


text:
Merge data frame on either of the 2 columns in R or python
I have two data frames. One data frame with one column and second data frame with two columns. I need to merge first data frame with either column of the second data frame and returns the values. After merging the values,the order of the first column should be same as the input. Please find the sample input below
first_column = c("ENSG00000165588","ENSG00000213551")
df1 = data.frame(first_column)
first_column = c("ENSG00000142192", "ENSG00000140575", "ENSG00000165588", "ENSG00000165588", "ENSG00000213551", "ENSG00000213551","ENSG00000197153")
second_column = c("ENSG00000165588", "ENSG00000165588", "ENSG00000186908", "ENSG00000135446", "ENSG00000273983", "ENSG00000274267","ENSG00000213551")
df2 = data.frame(first_column,second_column)

Please find the output below.
first_column = c("ENSG00000165588","ENSG00000165588","ENSG00000165588","ENSG00000165588","ENSG00000213551","ENSG00000213551","ENSG00000213551")
second_column = c("ENSG00000142192","ENSG00000140575","ENSG00000186908","ENSG00000135446","ENSG00000273983","ENSG00000274267","ENSG00000197153")
output = data.frame(first_column,second_column)

Just to confirm are you happy for a Python pandas solution? I recognise this is the implication of the title, but the code is in R

Yes, I am happy with python

This does not look like something I'd call a data frame merge. Are you able to give a bit more context around this problem (i.e. what you want to do with this data, rather than the intermediate problem you're trying to solve)? This helps us to wrap our heads around the problem and give better answers.

History This question was shifted from 'deleted' to 'closed' by me, which apparently was a mistake - apologies therein and then moved from 'closed' to 'opened'. I've now fully addressed the question. There was a misunderstanding: I originally provided a column merge solution when the OP wanted an index merge. Both solutions have now been provided. All code is vectorised and will be fast over a large data set.

I think this question is more about subsetting than merging. You can try this solution in R:
first_column = c("ENSG00000165588","ENSG00000213551")
df1 = data.frame(first_column)
first_column = c("ENSG00000142192", "ENSG00000140575", "ENSG00000165588", "ENSG00000165588",
                 "ENSG00000213551", "ENSG00000213551","ENSG00000197153")
second_column = c("ENSG00000165588", "ENSG00000165588", "ENSG00000186908", "ENSG00000135446",
                  "ENSG00000273983", "ENSG00000274267","ENSG00000213551")
df2 = data.frame(first_column,second_column)

output = rbind(df2[df2$first_column %in% df1$first_column,],
               df2[df2$second_column %in% df1$first_column,])

output <- transform(output, second_column = ifelse(output$second_column%in%df1$first_column,
                                                   first_column, second_column),
                            first_column = ifelse(output$second_column%in%df1$first_column,
                                                  second_column, first_column))
output <- output[ order(row.names(output)), ]

Edit The OP in hindsight is requesting an index merge. The first solution, presented here as Column merge solution was rejected by the OP.
From the comments below this post the order requested by the OP was:
# 0 ENSG00000165588 ENSG00000142192
# 1 ENSG00000165588 ENSG00000140575
# 2 ENSG00000165588 ENSG00000186908
# 3 ENSG00000165588 ENSG00000135446
# 4 ENSG00000213551 ENSG00000273983
# 5 ENSG00000213551 ENSG00000274267
# 6 ENSG00000213551 ENSG00000197153,

Column merge solution

The first operation is a concatenation and inner merge *.
The second operation are two outer merges, the second requiring a 'left-hand' merge *.

The code is here:
import pandas as pd

df1 = pd.DataFrame({'A':["ENSG00000165588","ENSG00000213551"]})
df2i = pd.DataFrame({'A':["ENSG00000142192", "ENSG00000140575", "ENSG00000165588",\
                          "ENSG00000165588", "ENSG00000213551", "ENSG00000213551", "ENSG00000197153"]})
df2ii = pd.DataFrame({'A':["ENSG00000165588", "ENSG00000165588", "ENSG00000186908", \
                            "ENSG00000135446", "ENSG00000273983", "ENSG00000274267","ENSG00000213551"]})

# Column 'A'
df3in = pd.concat([df2i, df2ii])
df4inner = pd.merge(df3in,df1, on='A', how='inner')

# Column 'A+'
df3out = pd.merge(df2i, df2ii, on='A', how='outer')
df4outer = pd.merge(df3in,df1, on='A', how='outer', indicator = True).query('_merge=="left_only"')
df5 = df4outer.drop('_merge', axis=1).reset_index(drop=True).rename(columns={'A':'A+'})
df6 = df4inner.join(df5)

print (df6)

Output
                 A               A+
0  ENSG00000165588  ENSG00000142192
1  ENSG00000165588  ENSG00000140575
2  ENSG00000165588  ENSG00000197153
3  ENSG00000165588  ENSG00000186908
4  ENSG00000213551  ENSG00000135446
5  ENSG00000213551  ENSG00000273983
6  ENSG00000213551  ENSG00000274267

This approach was rejected (comments) because of the position of ENSG00000197153 at row 2, Column A+ in relation to ENSG00000165588.

Index merge solution
import pandas as pd

def indexmerge(dfa, dfb, df1):
    df3out = dfa.reset_index().merge(dfb, how="left").set_index('index')
    df4outer = pd.merge(df3out,df1, on='A', how='outer', indicator = True).query('_merge=="left_only"')
    return df4outer

if __name__ == '__main__':

    df1 = pd.DataFrame({'A':["ENSG00000165588","ENSG00000213551"]})
    df2i = pd.DataFrame({'A':["ENSG00000142192", "ENSG00000140575", "ENSG00000165588",\
                              "ENSG00000165588", "ENSG00000213551", "ENSG00000213551", "ENSG00000197153"]})
    df2ii = pd.DataFrame({'A':["ENSG00000165588", "ENSG00000165588", "ENSG00000186908", \
                                "ENSG00000135446", "ENSG00000273983", "ENSG00000274267","ENSG00000213551"]})
    df3in = pd.concat([df2i, df2ii])
    df4inner = pd.merge(df3in,df1, on='A', how='inner')

    df4outer = indexmerge(df2i,df2ii,df1)
    df4outer2 = indexmerge(df2ii,df2i,df1)
    df4outer2.rename(columns={'A':'_A'},inplace=True)
    df6 = pd.concat([df4outer,df4outer2],axis=1)
    df6['A+'] = df6['A'].astype(str).replace('nan','') + df6['_A'].astype(str).replace('nan','')
    df7 = df6.sort_index().reset_index()
    df8 = df4inner.join(df7['A+'])
    print (df8)

Output
                 A               A+
0  ENSG00000165588  ENSG00000142192
1  ENSG00000165588  ENSG00000140575
2  ENSG00000165588  ENSG00000186908
3  ENSG00000165588  ENSG00000135446
4  ENSG00000213551  ENSG00000273983
5  ENSG00000213551  ENSG00000274267
6  ENSG00000213551  ENSG00000197153

This solution is identical to the order requested in the comments and is now also presented at the beginning of this post. The solution used here pivots around line 1 of indexmerge() and thus is an index merge operation. The remaining code is common pandas manipulation supporting this operation.
Note the new position of ENSG00000197153 row 6 column A+. The index merge retains ENSG00000197153 at the base of the column and therefore in association, i.e. same row, as ENSG00000213551.

* technical pandas terminology

Can you please check the output. The output must be in the dataframe. The order may be different but column A and column B entries should match. "ENSG00000213551" match in one case is wrong. If possible check my sample output above.

0 ENSG00000165588 ENSG00000142192
2 ENSG00000165588 ENSG00000140575
3 ENSG00000165588 ENSG00000186908
4 ENSG00000165588 ENSG00000135446
5 ENSG00000213551 ENSG00000273983
6 ENSG00000213551 ENSG00000274267
7 ENSG00000213551 ENSG00000197153, this is the output which i need.  For eg. ENSG00000213551 should interact with ENSG00000273983, ENSG00000274267, ENSG00000197153, but your program output shows ENSG00000213551 interact with ENSG00000135446, ENSG00000274267,  ENSG00000197153.

I have added the sample output and my question is correct.

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


text:
Sanger sequencing annotation error
I am a student in a Cancer lab. Working with sanger is new to me. While analyzing a report we found an insertion that has not been reported in any databases so far, we were working on checking if the insertion is pathogenic or not but when we reanalyzed the sanger annotation through variant reporter of the insertion it was giving it one position ahead and we analyzed that with the VEP as well that is were the confusion started.
The sequencer showed annotation of an insertion not between the bases but a step/position next to it:
c.5532-5533insGAA is the insertion position came after position 32 and before 33 but the sequencer annotates it like c.5533-5534insGAA.
We are working on human genome hg38.
Should I consider it a machine error or should I take that position for the result? Or is there a nomenclature error I am making here?
And
the annotation is in this form [c.33-34insAUU]+[=] and I don't know what the +[=] represents here.

How do you know it should be c.32-33insAUU and not c.33-34insAUU? Could this be an issue of 0-based vs 1-based position counting?

I haven't heard of the 0-based and 1-based system. I read the linked post and it would make sense if that is the case. The insertion in that case would be 0-based, I assume. But it would be hard to say if that really is the case since I don't about this system.
And it would possibly affect the result analysis, right? bcz the 1st annotation in my data tells its frameshift insertion while the 2nd annotation shows protein termination even when there is no stop codon there, according to Ensembl VEP results.
(the annotation in the question are examples)

Please [edit] your question and add more detail. You said you were performing Sanger sequencing, but then mention VEP. VEP requires specific variants, so can only be used after sequencing. Please explain what you are really doing so we can understand. Basically, who is telling you c.32-33insAUU and who is telling you c.33-34insAUU? What input data is given in each case?

Please EDIT your question and explain what you have. Sanger sequencing doesn't give variants, it gives sequences. So presumably you are comparing to some reference genome and inferring the HGVS notation, right? Tell us. Then you are using VEP? OK that needs to be in the question. If possible, you should also give us the actual, full HGVS notation so we can look it up and see what the reference genome has in that position. I am guessing you are working with human? Which genome? hg19? hg38? Did VEP use the same one as you did originally?

Actually, I dont have full information on that right now. I'll have to look it up and be more specific in the details of the questions and let me edit it. Thank you

Try looking up the variant on a tool like https://varsome.com (disclaimer: I work for the company that builds that tool, but it's free) where you can see the underlying genome sequence, and can easily switch between hg19 and hg38. VarSome will also tell you if there are multiple, equivalent indels (various variants that all cause the same genotype, think of a deletion of A from AAA: no matter which A you delete, you end up with AA).

Thank you for your guidance, I'll redo my post once I gather more specific information.

Just to give my thoughts:

The missing information is the actual sequence

Is this a conflict in the change in position of the insert?
Could it simply be a consequence of down stream numbering, either a logical error (coding stuff), or lab artefact.

My bet is its point 2, but really we don't know except for the OP because we don't know the sequence.
Sources of error, Point 2

Sanger sequencing was traditionally known for "frame-shift" errors if there were mono-mers, e.g. AAAA -> AAAAA, it can mess up over the number of bases in a mono-mer.
PCR amplification can also cause slippage against monomers, same thing as the above point - but it's real.
1 vs 0 ... this is because in a Pythonic list or any C/Perl-based array numbering starts at 0 (not 1) - hence @terdon's comments.

If upstream mono-mers are present that might be the source of the confusion here and given the amount of sequence data thats very likely.
What I'd do:

Just Blast the new sequence and see what emerges: if it is genuinely novel - they'll be zero 100% identity hits.
VCFs ... I think this is VarSome's stuff (forget).

If there are Blast hits it will identify whether there is any upstream lab artefacts.
Conclusion
If this was any other situation in human genetics, or even mammalian genetics I'd be excited. However, in this case - I think its a routine result. Firstly, the chance of finding a new "driver" mutation are slim. Secondly, I've looked at this chromosome in a cancer patient (only once) and it was nuts (and thats an under statement) - I've a fairly good memory for numerical patterns, and I'd never seen a pathogenic species behaving like that. I think this insert is just a consequence of the entire mutational chaos that are cancer cells.
Just to make it clear: I don't do cancer but if you want an authoritative assessment and think you've something major I'd contact (or your PI) the Broad Institute via Geraldine A. Van der Auwera who also has an account on this site via @Geraldine_VdAuwera. Sanger sequencing is not what they do, but they'd know at the drop of a hat what the significance and if they don't respond - thats your answer ;-). Alternatively, VarSome would know, albeit they're commercial enterprise offering freeware.

@terdon
I have the reads generated by sanger. If there is a tool I can use to get the insertion coordinates to confirm? The sanger generated  results are doubtful because according to them the protein terminates due to insertion between a codon that would indicate severe condion which is not suggested in the patient's history.

The insertion is GAA so is it possible that the sequence shifted the results one base ahead? I have a picture of graph if you want to have alook

One thing you can try, @sanaamir, (and, again, I am pushing my own company's product here; I'm sorry, but those are the tools I am most familiar with) simply pasting the read into VarSome. That will align it against hg19 or hg38 (your choice), and identify any variants found in that read. That will help you see the full notation of the variant.

Thank you I will do that. But do i put the files from sanger or just fasta format of nucleotides with the insertion? I did have a look at the website today. @terdon

Please come into the chat room linked to the question (https://chat.stackexchange.com/rooms/141809/discussion-between-sana-amir-and-terdon) so we don't spam M_.


text:
Add or Simulate Noise in Short-Read Paired End Data
For a project I'm working on, I need to figure out how to model noise that may be happening in real genomes due to alignment errors, contamination, etc. Specifically, short-read paired-end data either the fastqs or bams. The reason I want to do so is to distinguish this noise from signals of actual subclonal variants.
When I say noise, I mean things like discordant pairs with abnormal orientation happening spuriously across the genome. This noise might be happening due to alignment error or something else, I'm not entirely sure yet.
I've simulated data with SVEngine but the amount of noise (discordant pairs with 1 kb+ insert sizes) that don't overlap with simulated variants is minimal. A total of 11 pairs.
Is there a tool or method I can use to add this noise, specifically discordant pair noise? I don't want to add subclonal variants with a low VAF as a replacement for noise, rather add something truly spurious.
Thanks in advance!

is what you want a true simulation approach or a sampling approach? my (unfounded) impression is that a lot of the idiopathic technical error that you describe is something that is not captured by error models. however, by drawing directly from real read data you will likely capture this error profile, though at the cost of including whatever real variation there is. of course, you will want to first confirm that no simulator accurately captures the dynamics that you want. possibly see this benchmarking study.

Possible duplicate of this question


text:
How to calculate the Maximum Ladder Distance (MLD) of ssRNA from Dot-Bracket notation?
The ladder distance (LD) between residue i and residue j is the minimum number of base pairs (created due to secondary structure of the ssRNA) one has to cross when going from i to j.
The Maximum Ladder Distance of a ssRNA molecule is the maximum ladder distance between any two of the residues of this molecule[*].
Given the secondary structure of a ssRNA molecule in dot-bracket notation (e.g. by RNAfold), how can I calculate the Maximum Ladder Distance of this molecule?
Below is an illustration from [*]. The MLD is 11, since this is the number of basepairs on the yellow line.
Edit: The dot-bracket notation for this molecule would be
..(((((.......))...((.((((((.......))))))..)).)).

*: Yoffe AM, Prinsen P, Gopal A, Knobler CM, Gelbart WM, et al. (2008)
Predicting the sizes of large RNA molecules. Proc Natl Acad Sci U S A 105:
16153–16158

Please [edit] your question and post the dot-bracket representation of the above structure. Use the code formatting button ({}) to format it in a monospace font for clarity.

There are different implementations for calcuating the $\rm MLD$. One possible approach is to consider $\rm MLD$ as the heighest value of all possible mountain plots, generated by shifting the beginning of your sequence. The highest value of a mountain plot tells you the maximal ladder distance that can be reached from the first nucleotide. But the longest path might not begin at the first nucleotide, and that is why you have to go through all possible shifts of the sequence.
Here is an example of this implementation in Python. The algorithm is as following:

Convert dot-bracket notation to a list of base pairs (dot2pairs).
Calculate the first mountain plot.
Calculate required shifts. This is not strictly neccessary because you can swipe through all possible shifts = list(range(0, length)), but it gives better performance.
Iterate through all shifts, calculate the MLD at each step, and return the maximal one.

import numpy as np
from collections import deque

def dot2pairs(dot):
    stack, pairs = deque(), []

    for i, n in enumerate(dot):
        if n == '(':
            stack.append(i)
        elif n == ')':
            pairs.append((stack.pop(), i))

    return np.array(pairs)

def dot2MLD(dot):
    # Convert to a list of base pairs
    pairs0 = dot2pairs(dot)

    if len(pairs0) == 0:
        return 0

    MLD = 0
    shifts = [0]
    length = len(dot)
    seq = np.zeros(length, dtype=int)
    curr_mount = None

    while shifts:
        # Shift the sequence
        pairs = (pairs0 - shifts.pop()) % length

        # Find positions of opening and closing nucleotides
        opens, closes = np.min(pairs, axis=1), np.max(pairs, axis=1)

        seq.fill(0)
        seq[opens] = 1
        seq[closes] = -1

        mountain = np.cumsum(seq)

        MLD = max(MLD, int(np.max(mountain)))

        # In first iteration, calculate all required shifts
        if curr_mount is None:
            curr_mount = mountain[0]
            for i, m in enumerate(mountain):
                if m == curr_mount:
                    shifts.append(i)
                else:
                    curr_mount = m

    return MLD

Let's test it on your structure:
dot2MLD("..(((((.......)))...((.((((((.......))))))..)).)).")

11

And on a larger structure:
dot2MLD(s)

1141

It took ~ 5 seconds for this 17481 nt long structure (and it takes ~ 20 seconds if you don't precalculate requires shifts).


text:
Do you need galaxy.yml when you start creating galaxy?
First question: I found that there is no galaxy.yml in the config directory of the running version of docker galaxy, only galaxy.yml.sample file and galaxy.yml.interactivatetools file, so do I need to create it, or do I need to create a new one using the cat command? galaxy.yml file?
The second question is about these initial notes, where can I go deeper?

What have you already tried? Does Galaxy work for you if you run it without a .yml file? Have you looked at the galaxy documentation on the Galaxy website?


text:
Is there available REGA or COMET versions for downloading & using on a local machine?
I am currently working on verifying the subgenotyping of our sequences and a subsample from the Los Alamos NL database. My supervisor has recommended using REGA & COMET for this purpose. While our dataset, which consists of slightly > 1,000 sequences, can be easily subgenotyped using the online versions of these tools, the LANL dataset is significantly larger, prompting me to consider local installations. Unfortunately, however, I have been unable to find the source repositories for REGA & COMET.
I did come across a link to the entire REGA DB but I am unsure if it contains the subgenotyping tool or if the last can be installed separately. I am also interested in exploring other locally installable tools of comparable quality.
In my search, I found jpHMM, which is available for download and free use. However, it seems to be less frequently used in publications compared to the aforementioned ones.
I would appreciate any insights on the quality of jpHMM or any other recommended tools for local installation (even as packages, e.g. for Python).
If you have experience with these subgenotyping tools or any relevant advice, please share your thoughts! Your help is greatly appreciated!
Best regards,


text:
Automated method to figure out whether a recognition site is conserved or not using the results of simulated restriction digestion

Restrict from EMBOSS Suit version 6.3.1 with the following
parameters: snucleotide1, sitelen = 4, rformat = table, enzymes =
enzymes.txt.
From the 4379 enzymes present in REBASE, we selected the 650 restriction enzymes that were commercially available, since this assay is meant to be used in any laboratory. From the 650 enzymes, 152 digest all P. salmonis sequences and only 65 recognized conserved restriction sites in the complete set of sequences, generating the same/similar restriction pattern (same number of bands and similar sizes).

https://www.frontiersin.org/articles/10.3389/fmicb.2016.00643/full
If my reading is correct, these researchers figured out the conservation of the restriction sites from the resulting virtual patterns.
The coordinates of a cut can be different in two related sequences even if the restriction site is conserved. This is because base insertions/deletions could occur before/after the restriction site.
How do i figure out if a recognition sequence is conserved or not, if the coordinates of the cut are not the same in the sequences? How do i automate this?
We can consider this as an example:
I have two sequences of the 16s rRNA gene from two strains of a bacterial species. When i run EMBOSS's restrict by the following command:

restrict -sitelen 4 -enzymes HaeIII -snucleotide_sequence 1 -sequence
sequences.fasta -rformat excel -fragments true -outfile patterns.xls

Sequences link
The result is:


text:
Filter VCF file for most differentiated SNPs (Fst) in VCFTools
I have a filtered VCF file that I want to subset for the 500 most differentiated SNPs based on estimated Fst values (calculated in vcftools). I am not completely sure how to approach this. I have calculated the fst values across all populations (10 in total), but am not sure how to specify that I only want 500 SNPs with the highest fst value.
I used VCFTools, but would BCFTools be more appropriate?

Please [edit] your question and give us an example we can use, and the output you would want from that example. We cannot help you filter data we cannot see. Show us a valid VCF file (that means all headers plus a few variant lines) and then tell us which variants you want to keep from it.


text:
overrepresentation test, between transcriptome and candidates sequences obtained from the transcriptome
For an analysis of my data, I have a transcriptome and a list of sequences obtained from the transcriptome. I would like to perform a functional enrichment analysis. I have annotated both sets of data using eggnog mapper. Currently, I want to perform a test between the two functional annotations, specifically COGs (Clusters of Orthologous Groups). I have tried using the R code https://yulab-smu.top/biomedical-knowledge-mining-book/enrichment-overview.html#gsea-algorithm
with clusterProfiler, but it seems that it may not work. With which tools or code can I perform this test, please?
exemple somme of my data :


text:
How to find the correlation between two genes and their separate overall survival outcome?
i am struggling to find the correlation between two isoforms and their survival outcome. i do have a dataset that is a cohort for thousands of genes and am only interesting in two which they should be isofroms to each other. However, i thought the ratio between these two might be a prognostic marker. i need to find a way where i can see if this is possible by interpreting the clinical outcome which is the overall survival. many thanks guys.
Crossposted: https://www.biostars.org/p/9566172/

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.


text:
adding multiple files in folder and zipping in loop
Through fasterq-dump, I am generating 2 files in 1 run of loop. I want to club these two files (A_1.fastq & A_2.fastq) in one folder and zip the it.
cat SRR_Acc_List.txt | while read i ; do  fasterq-dump ${i} --split-files -v --skip-technical  -e 15 ; for j in 2 ; do   gzip ${i}_j.fastq  ; done; done
cat SRR_Acc_List.txt | while read i ; do  fasterq-dump ${i} --split-files -v --skip-technical  -e 15 ; for(j = 1;j<2;j++) ; do   gzip ${i}_j.fastq  ; done; done

cat SRR_Acc_List.txt | while read i ; do  fasterq-dump ${i} --split-files -v --skip-technical  -e 15 ; for j in 1 2 ; do   gzip ${i}_${j}.fastq  ; done; done

cat SRR_Acc_List.txt | while read i ; do  fasterq-dump ${i} --split-files -v --skip-technical  -e 15 ; for j in 2 ; do   gzip ${i}_${j}.fastq  ; done; done
cat SRR_Acc_List.txt | while read i ; do  fasterq-dump ${i} --split-files -v --skip-technical  -e 15 | mkdir {i} | for j in i; do   gzip ${i}  ; done; done

cat SRR_Acc_List.txt | while read i ; do  fasterq-dump ${i} --split-files -v --skip-technical  -e 15 | mkdir ${i} ; for j in i; do   gzip ${i}  ; done; done
cat SRR_Acc_List.txt | while read i ; do  fasterq-dump ${i} --split-files -v --skip-technical  -e 15 ; mkdir ${i} ; mv ${i}_*.fastq "${i}"; for j in i ; do gzip ${i} ; done;done

cat SRR_Acc_List.txt | while read i ; do  fasterq-dump ${i} --split-files -v --skip-technical  -e 15 ; mkdir ${i} ; mv ${i}_*.fastq "${i}";done; for j in i ; do zip ${i} ; done

cat SRR_Acc_List.txt | while read i ; do  fasterq-dump ${i} --split-files -v --skip-technical  -e 15 ; mkdir ${i} ; mv ${i}_*.fastq "${i}" ;  zip ${i} ; done

cat SRR_Acc_List.txt | while read i ; do  fasterq-dump ${i} --split-files -v --skip-technical  -e 15 ; mkdir ${i} ; mv ${i}_*.fastq "${i}" ; tar -cv ${i} | gzip > ${i}.tar.gz ; done

None of them is working.

Try your code on one set of FQs and use a loop only after that set gets processed successfully. Also, why pipe tar to gzip when you can simply do tar -czf archive.tar.gz files_to_archive?


text:
Compare in Python / PyRanges two genomic intervals to determine whether they're separate, partially overlapping, or contained
As part of a larger analysis, I need to evaluate pairs of genomic regions (regionsA and regionsB) and perform different calculations for each pair depending on whether:

They're completely separate (no overlap)
regionsB partially overlaps regionsA (start or end of regionsB does not overlap; both are handled separately)
regionsB is fully contained in regionsA

At this time, regionsA and regionsB are represented as Pandas DataFrames. Likewise, I'm doing all possible combinations:
if (region_b.start < region_a.start and
    region_b.end < region_b.start):
[...]
elif (region_b.start < region_a.start and
     (region_b.end > region_a.start
     and region_b.end < region_a.end )):

Given that there are many conditions I check for, this makes the code long, a little hard to read, and too verbose. Is there anything I can use (without making the code harder to read) to make these checks simpler, something like PyRanges (or PyRanges itself, although I don't think there's API for it)?
If not, I'll just live with what I have.


text:
Error while opening xmfa generated by progressiveMauve
I used commandline progressiveMauve to align two genome .gbk files. The alignment was successful. But when I try to visualize the final .xmfa file with Mauve, I get the following error.
Exception in thread "Thread-1" java.lang.RuntimeException: Similarity index is too large.
at org.gel.mauve.histogram.ZoomHistogram.allocateIndex(Unknown Source)
at org.gel.mauve.histogram.ZoomHistogram.<init>(Unknown Source)
at org.gel.mauve.SimilarityIndex.<init>(Unknown Source)
at org.gel.mauve.XmfaViewerModel.init(Unknown Source)
at org.gel.mauve.XmfaViewerModel.<init>(Unknown Source)
at org.gel.mauve.ModelBuilder.buildModel(Unknown Source)
at org.gel.mauve.gui.FrameLoader.loadFile(Unknown Source)
at org.gel.mauve.gui.FrameLoader.run(Unknown Source)
at java.lang.Thread.run(Thread.java:745)

How can I solve this?
Is there another way to visualize the xmfa file.
My aim is to visualize an inversion in one of the genomes.


text:
Using PANDA to infer gene regulatory networks without ppi
I'm using this tutorial (https://github.com/netZoo/netZooM/blob/master/tutorials/panda/panda.pdf) to use PANDA to infer gene regulatory networks, but without ppi values (which I've been told can be done easily using MATLAB). I replicated the codes given in the tutorial but I omitted all the portions with ppi. After doing so, I'm getting a lot of warnings and errors that makes me think I've made syntax errors in the process of taking out the ppi portions. Here's what I have:
addpath(genpath('path/to/netZooM'))

load("myMat.mat")
myTable = struct2table(myTable);
addpath(genpath('path/to/netZooM'));
fileparts(which('netZooM'));

exp_file = 'tests/test_data/expression.txt';
motif_file = 'tests/test_data/motifTest.txt';
panda_out = '';
lib_path = 'C:/Users/ual-laptop/Documents/bioinformatics.m/netZooM';
save_temp = '';

alpha = 0.1;
save_pairs = 0;

expressionData = readTable('expression.txt');
size(expressionData)

[TF, gene, weight] = textread(motif_file, '%s%s%f');
TFNames = unique(TF);

disp('Reading in expression data!');
fid = fopen(exp file, 'r');
headings = fget1(fid)
n = length(regexp(headings, '\t'));
frewind(fid);
Exp = textscan(fid, ['%s', repmat('%f', 1, n)], 'delimiter', '\t');
fclose(fid);
GeneNames = Exp{1};
Exp = cat(2, Exp{2:end});
[NumGenes, NumConditions] = size(Exp);
fprintf('%d genes and %d conditions!\n', NumGenes, NumConditions);
Exp = Exp';

disp('Computing coexpression network:');
GeneCoReg = Coexpression(Exp);

disp('Reading in motif data!');
[TF, gene, weight] = textread(motif_file, '%s%s%f');
TFNames = unique(TF);
NumTFs = length(TFNames);
[~,i] = ismember(TF, TFNames);
[~,j] = ismember(gene, GeneNames);
RegNet = zeros(NumTFs, NumGenes);
RegNet(sub2ind([NumTFs, NumGenes], i, j));
fprintf('%d TFs and %d edges!\n', NumTFs, length(weight));

disp('Normalizing Networks:');
RegNet = NormalizeNetwork(RegNet);
GeneCoReg = NormalizeNetwork(GeneCoReg);
TFCoop = NormalizeNetwork(TFCoop);

disp('Running PANDA algorithm:');
AgNet = PANDA(RegNet, GeneCoReg, TFCoop, alpha);

AgNet = panda_run(lib_path, exp_file, motif_file, panda_out, save_temp, alpha);

sum(AgNet, 2)
sum(Agnet, 1)

Are there any glaring syntax errors in my code? I'm thinking my errors could also stem from maybe cloning the github repository into the wrong place (which is a separate topic) but it's definitely possible that there are syntax errors since this is my first time doing this kind of stuff.


text:
Adding BAF and Log2 ratio to VCF file
I have access to pre-processed VCF files from a genetic dataset, and am trying to use the MoChA tool to calculate mosaic chromosomal alterations in this dataset. However, MoChA requires BAF and log2ratio values for each of the variants, which are missing from the VCF files given, and the original genome array CEL files are unavailable so I cannot redo the variant calling and include these statistics
Instead, they have provided us with nexus format CNV files (from the Axiom CNV summary tools) for each of the samples with a host of log2 ratio and BAF files across the different Axiom probes.
I am just a bit confused about how we can use these files to import the BAF and log2ratio values into the VCF for MoChA. Any suggestions?
A couple sample rows of the CNV file are included here (I've removed the sample identifier and replaced with ---).
Sample  ProbeSetName    Chromosome      Position        Log2Ratio       BAF
---     AFFX-KIT-000001 9       101258881       -0.440535940556599      0.482727669963674
---     AFFX-KIT-000002 4       164934874       -0.106227164241716      0.057035393165642

Welcome to the site. The manual from GitHub https://github.com/freeseek/mocha appears to describe this in detail both for BAF and log2 (which the manual calls LRR for Log R Ratio). I would recommend reading the manual because it is extensive and looks very clear.


text:
Plot error and blank tab of plot
I am trying to plot by using these commands , but this is the error I am getting, and I tried multiple times, even by making the window size larger, the [lot tab is empty/blank too.
plot1 <- VariableFeaturePlot(nsclc.seurat.obj)
LabelPoints(plot = plot1, points = top10, repel = TRUE)***
*When using repel, set xnudge and ynudge to 0 for optimal results
Error in grid.Call(C_convert, x, as.integer(whatfrom), as.integer(whatto),  :
  Viewport has zero dimension(s)
In addition: Warning message:
Transformation introduced infinite values in continuous x-axis*

Welcome @Samia. Could you kindly clarify your study? I get the impression this is a Seurat analysis and I assume you are plotting the PCA. Any further information on the study would help understand the error output. Basically that the PCA is empty "zero dimension(s)" - not good. The error is likely further upstream.


text:
How does BioJava retrieve bond information from a PDB file?
How can I compute energy from a protein PDB file?
I asked the above question.
Then I thought, Why don't I check how BioJava obtains bond information from a protein PDB file?
Surprisingly, BioJava's tutorial doesn't have any discussion on Bonds.
I checked BioJava's GitHub repo, but I failed to find any such information.
How does BioJava retrieve bond information from a PDB file?


text:
Calculating p-value and adjusted p-value for pre-normalized microarray data with fold change precalculated
I have a microarray dataset with two mutants dataset that has already been normalised, and the fold change values for each gene in each mutant versus the wild type have been calculated. I'm interested in determining the statistical significance of these fold changes.
What I've tried: I initially considered using limma, but I wasn't sure if it was appropriate for pre-calculated fold change values.
What I was expecting: I was hoping to obtain a list of genes with their corresponding p-values and adjusted p-values, indicating their significance. Ideally, this would help me identify genes that are significantly upregulated or downregulated in the mutant compared to the wild type.
Given that I only have the fold change values for two biological replicates (and no raw intensity or count data), what would be the best approach to achieve this? Any guidance, recommended tools, or R packages would be greatly appreciated!
Thank you in advance!

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.


text:
Problem with edit and align viral whole genomes
I have a database of whole genome of African swine fever virus (about 190kb/genome) downloaded from GenBank. But these sequences do not start at the same place on the genome. When I align the sequence to build a phylogenetic tree. The sequences at the 3' and 5' ends have many gaps and do not match. The length of the entire data set is over 220kb due to many gaps.
Kindly suggest me a solution to the problem.

Welcome! Can you post an image of your alignment?


text:
Need to make circos plot from Enrichment data
I have a list of genes and their associated pathways. The data was generated using Partek. The data looks like this -

Gene_in_list
Gene_set
Description
Enrichment score
P.value

Gene1
hsa
Disease C
69
0.005

Gene1
hsa
Disease A
45
0.05

Gene2
hsa
Disease B
23
0.007

Gene2
hsa
Disease Y
40
0.009

Gene2
hsa
Disease C
69
0.0001

I want to 'split open' the Description column (Which I think means going long to wide in tidyverse terminology?), and transform it to the below table -

Gene_in_list
Disease A
Disease B
Disease C
Disease Y

Gene 1
1
0
1
0

Gene 2
0
1
1
1

So that I can use the online circos plot to essentially show, which genes are involved in the pathways.
I was also debating if Enrichment score / p.value could be used as the thickness of connections in circos, but I think it might look too thin / uneven.
Any help appreciated, thanks!


text:
Unable to access error logs of failed cromwell execution (WDL) in local machine running in ubuntu environment
I was trying to run bases2fastq software downloaded from the element bioscience platform. I was running this software through WDL workflow named list_fastq. For some reason, maybe syntax, my cromwell execution failed. However, I haven't been able to read the logs as they are empty. Below is the content of the 'stderr.background' file.
mkfifo: cannot create fifo '/mnt/d/element-bioscience/cromwell-executions/list_fastq/5e9a780d-c010-41f7-bba0-bcb42c870891/call-list_samples/tmp.05c20794/out.78994': Operation not supported
mkfifo: cannot create fifo '/mnt/d/element-bioscience/cromwell-executions/list_fastq/5e9a780d-c010-41f7-bba0-bcb42c870891/call-list_samples/tmp.05c20794/err.78994': Operation not supported
/mnt/d/element-bioscience/cromwell-executions/list_fastq/5e9a780d-c010-41f7-bba0-bcb42c870891/call-list_samples/execution/script: line 17: /mnt/d/element-bioscience/cromwell-executions/list_fastq/5e9a780d-c010-41f7-bba0-bcb42c870891/call-list_samples/tmp.05c20794/out.78994: No such file or directory
/mnt/d/element-bioscience/cromwell-executions/list_fastq/5e9a780d-c010-41f7-bba0-bcb42c870891/call-list_samples/execution/script: line 18: /mnt/d/element-bioscience/cromwell-executions/list_fastq/5e9a780d-c010-41f7-bba0-bcb42c870891/call-list_samples/tmp.05c20794/err.78994: No such file or directory

Please guide me to the steps on being able to prevent cromwell from deleting the contents of log files, if that's what's happening.


text:
PyRanges: drop specific record from an object
Given a specific PyRanges object, such as
import pyranges as pr
gr = pr.from_dict({"Chromosome": [1, 1, 1, 1], "Start": [1, 25000, 30000, 40000],
                   "End": [24000, 29500, 37000, 45000], "ratio": [0.5, 0.3, 0.5, 0.4],
                   "ID": (0,1,2,3), "Strand": "+"})

gr

+--------------+-----------+-----------+-------------+-----------+--------------+
|   Chromosome |     Start |       End |       ratio |        ID | Strand       |
|   (category) |   (int64) |   (int64) |   (float64) |   (int64) | (category)   |
|--------------+-----------+-----------+-------------+-----------+--------------|
|            1 |         1 |     24000 |         0.5 |         0 | +            |
|            1 |     25000 |     29500 |         0.3 |         1 | +            |
|            1 |     30000 |     37000 |         0.5 |         2 | +            |
|            1 |     40000 |     45000 |         0.4 |         3 | +            |
+--------------+-----------+-----------+-------------+-----------+--------------+

I want to get rid of one interval here (not necessarily the first or the last: for this exercise, for example the third one, without relying on the presence of the ID column.  So I would need to drop the record with Chromosome 1, Start 30000 and End 37000. The use case is that I need to replace it with a different record.
Replacement of the record is somewhat expensive but not a problem (pr.concat followed by a sort to put everything in order), but I can't figure, after reading docstrings and searching around, how to remove a record.
I can only think of making the record I want to drop a second PyRanges object then using overlap() with invert=True to remove the record. Is there a more elegant and less convoluted way to do this?


text:
how to use CIBERSORT
I have seen in several articles, the use of CIBERSORT to calculate the proportion of immune system cells using gene expression data, I would like to try it but I have not found how to install it for R, I have tried to register on the website but it doesn't allow me to enter after introduing my e-mail and password.

Could you clarify the issue ?


text:
What versions of pip and python are appropriate for nsdpy?
I am currently trying to download all the requirements to make a COInr database, following this workflow; https://mkcoinr.readthedocs.io/en/latest/content/overview.html.
It requires a couple different third party programs to be able to run, one of them being nsdpy, which it says all it requires is cloning the repo from Github and then you're supposed to simply be able to use 'pip' to install nsdpy.
However I continued to receive this error;
""DEPRECATION: Python 2.7 will reach the end of its life on January 1st, 2020. Please upgrade your Python as Python 2.7 won't be maintained after that date. A future version of pip will drop support for Python 2.7. More details about Python 2 support in pip, can be found at https://pip.pypa.io/en/latest/development/release-process/#python-2-support
ERROR: Could not find a version that satisfies the requirement nsdpy (from versions: none)
ERROR: No matching distribution found for nsdpy""
Which makes it seem that my issue is with the version of our Python (I am using a linux system on a Mac). The version of Python that is says we should be using is 2.7, the current version I have downloaded is 3.7.7. But then I was reading through the nsdpy docs (https://nsdpy.readthedocs.io/en/latest/main.html#overview) which says that it requires a Python version of 3.8+. So I am a little confused, which version of Python to I truly need to be able to install nsdpy?
Is there a more up to date version of nsdpy?
If anybody else has had this problem and/or has a solution, please let me know!

To add comment that was lost in migration: I strongly suspect that the version of pip is at fault, as it seems like the issue is that python2.7 is being run, not that it is required.


text:
Estimating rate of evolution from a phylogeny
Any suggestions for the best way to estimate a metric for evolutionary rate of viral sequences of the same species?
I was using terminal branch length from a phylogeny and would quite like to try sub rates but i am really new to this

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.


text:
Making a matrix when samples that have missing data
I am trying to make a dataset for methylation analysis.
I called methylation using Methyldackel and there are various positions without methylation calls. This creates a problem because when I try to use for example, DMCHMM or MethylSig, to estimate differentially methylated regions, the initial step prompts an error:
Building BSData object.
Warning message:
In .merge_two_Seqinfo_objects(x, y) :
  Each of the 2 combined objects has sequence levels not in the other:
  - in 'x': NW_022883351.1, NW_022883409.1
  - in 'y': NW_022883010.1
  Make sure to always combine/compare objects based on the same reference
  genome (use suppressWarnings() to suppress this warning).

I used the same reference for all my samples, but when it tries to combine the files, it complains because some samples don't have data (missing) in some positions.
Is there a way to use R or Linux to fix this, and check those lines positions with missing data and add those (lines) positions in the files?
The files look like this:
chrBase chr     base    strand  coverage        freqC   freqT
NW_022882922.1.4856     NW_022882922.1  4856    F       2         0.00  100.00
NW_022882922.1.4866     NW_022882922.1  4866    F       2         0.00  100.00
NW_022882922.1.4889     NW_022882922.1  4889    F       4         0.00  100.00
NW_022882922.1.4892     NW_022882922.1  4892    F       4         0.00  100.00

so that all the files have "information" in all the positions?
Thank you
Juan Pablo


text:
How can I use QIAGEN Ingenuity Pathway Analysis (IPA) on proteomics data associated with an outcome?
I am trying to do pathway analysis on a set of proteins associated with several outcomes (one is continuous, the others four are binary). I have run independent associations between the proteomics data and the outcomes in question and have ascertained betas, hazard ratios, p values, and associated multiple-testing correction p-values. Out of the set of 3000 proteins for which we have proteomics data, we have found significant FDR corrected associations with our outcomes in ~1800. I am trying to do pathway analysis with these proteins.
I am trying several pathway analysis approaches, including fgsea, online GSEA through GenePattern, and Ingenuity Pathway Analysis. For the GSEA approaches, I have been able to generate the ranked list and provide the ranked list to the software. However, it is my understanding that IPA is a form of overrepresentation analysis and therefore I should not provide a ranked list. However, when uploading my dataset to the IPA application, it asks for expression change values (my understanding is that all these approaches were originally designed for microarray-based  differential expression data).
Since we are looking at the association between proteomics and several outcomes, expression change values are not relevant to these data. However, the directionality of the effect is relevant when considering pathways, and since we have the betas for the association it seems to me that I should include it. I looked through the IPA documentation and all of it seems specific to differential expression analysis. Can anyone point me in the correct direction to perform IPA on the association between my outcomes and proteomics data? Any help would be appreciated.

What have you already asked Qiagen about this? IPA is a very expensive product; Qiagen should at least provide good support around the bioinformatics side of things because bioinformatics is such a core part of their product.


text:
How do I get the gene annotations as a text file from a genbank file?
Also posted on biostars
I have a genbank file and I need the following details about the gene annotations in a text file
# ID    name    start   stop    frame   product description function    cogcats status_region   status_function type    region_number   function_number

How do I go about extracting them? Please let me know if more information is needed from me! Thank you

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.

Please [edit] your question and add i) an example input genbank file and ii) the exact output you would want from that example. Also tell us what operating system you are using so we know what tools you have available. Will this file describe a transcript? A gene? A chromosome? Most of the fields you mention aren't standard, so we really need the example to understand.

The fact that Genbank files are also text files, I believe you meant that you prefer the same data in the "table format". In that sense why not fetching your annotations of interest in GFF/GTF format? And if that is not available for your genome of interest, there are tools that convert Genbank to GFF. Seeing the column names in your example, I think you can go for the GFF3 format and "expand" the last column (attritubes) to extract the semicolon-separated info within that column.

Crossposted here: https://www.biostars.org/p/9591298/


text:
Dante Labs VCF: UPDATE: SOLVED
Has anyone tried to analyse the VCF files provided by Dante Labs? (Human WGS, Illumina)
I'm looking into .raw.vcf and .filtered.vcf.
For reference sequence, all I can see is 2 contigs:
GL000225.1 and GL000192.1.
All variants belong to one of these.
There's no conventional CHROM column.
I don't know how those coordinates map to chromosomes in GRCh37 or GRCh38 assemblies.

GL000192 corresponds to chr1_gl000192_random
GL000225 corresponds to chrUn_gl000225

Moreover, there's no VCF header with metainfo, all I can see is 1 top row
:11,27:0.711:38:4,11:7,16:25:64,0,24:2.8811e+01,1.4126e-02,2.7139e+01:0.00,34.77,37.77:6,5,12,15:7,4,18,9.

What's the use of those VCFs?
What variants do they contain?
What version of VCF it is and where can I get the standard metainfo?
Is there a map of those coordinates to any core assembly
chromosomes?
Am I right they contain not all genomic VCFs but only a small part
of the genome? (Size of GL000225.1 contig is 211,173 bp and size of
GL000192.1 is 547,496 bp, that obviously doesn't sum up to the whole
genome).

Many thanks for your replies!
UPD: I found a mistake, sorry. The file was incomplete, only tail was uploaded to Jupyterhub.
In local Jupyter, the file loaded completely, it actually has the  header and CHROM column.


text:
Generating synthetic RNA-Seq data for performance assessment of new algorithm
I want to generate synthetic miRNA data for the performance assessment of the newly developed algorithm designed for miRNA data analysis. There are several synthetic sequence simulators available in the published literature, such as wgsim, dwgsim, ART, inSillicoseq, mason,  NEAT, BEER, polyester, RSEM, SimSeq, etc. All of these tools are either for DNA or transcriptome profiles. Can anybody please suggest the following:

Can synthetic DNA sequence simulators (ART, wgsim, dwgsim) generate synthetic miRNA data?
Some tools, such as CAMPAREE, RSEM, and RNA-Seq simulator, generate RNA-seq data, but I think these tools generate mRNA data, and I am unsure whether I can use them for my miRNA data analysis.

Please suggest. I very much appreciate any help in this regard. Thanks


text:
Is there any tools to predict the potential RNA candidate that can bind to a protein?
I have a list of RNA binding proteins from HyDRa prediction. I have a long list of mRNA, I'm interested in.
For the list of RNA binding proteins, I want to find the potential RNA binding partners. Lets say each protein can bind to 5-20 RNAs, I want to find this list.
Now for each RNA-protein binding, I want to calculate the binding affinity or any equivalent measurement.
Is there any tools, I can use to do this? or any related insilico approach I can try?

Almost certainly I would look at published machine learning studies


text:
How to use CIRCOS plot pairwise syteny blocks for multi species
Hi, I want to use Circos to draw syntenic relationships among multiple species, like the Fig.1b in this article. How to define the ideograms for different species?
Thank you.
https://doi.org/10.1038/s41588-024-01683-

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.


text:
Any other (good) long read aligners apart kma?
I wonder if there are any other more efficient long read aligners apart kma? I have heard there is some hisat2 long read application but it's not really efficient. Any suggestions?

"Apps"? You mean aligners?

Yes, aligners. I shouldn't call them apps, I just meant some module that can be attached to a script, like bwa, minimap or that kma for ex.

"module that can be attached to a script" is not a thing. Modulefiles are a way of packaging binaries and related dependencies, so you're looking for an aligner like ATPoint says, not an app or a module. If you need a binary as a module, consult your HPC sysadmin.

Please [edit] your question and clarify your needs. If you need an aligner, then please say so (apps is meanigless). If you instead mean a variant caller, say that. Then, explain what more efficient means for you (faster? better resource utilization? less RAM specifically? something else?), and whether you are tied to a specific operating system. Perhaps also tell us what kind of sequencing you will be doing (targeted? WES? WGS?).