text:
Primer trimming-fasta files
My goal:
I want to trim off the primers (Forward : CGAGAAGACCCTRTGRAGCT, Reverse : GTTGGGGYGACCNYGG) from a fasta file with a lot of dna sequences allowing for some (e.g. 3) mismatches (identity).
Input sequences
Seq1:

CGAGAAGACCCTATGGAGCTTAAGGCGCCAGAACAGCTCACGTCAAACACCCCCGCATAAAGGGAATAAACCAAGTGGACCCTGCTCTAGTGTCTTTGGTTGGGGCGACCGCGGAACGT

Seq2:

ATGGCCCATCCCTCACAGCTAGGATTTCAAGATGCAGCTTCCCCAGTTATAGAAGAACTTCTCCACTTTCACGACCATGCCCTAATAATCGTTTTTTTAATTAGTACACTAGTA

Seq3:

CATAAGACGAAAAGACCCTATGGAGCTTTAGACGTCAGAGCAGCTCATGTAAAGCACCCCTAAACAAAGGAAAAAACCAAATGAAATCTGCCCTAATGTCTTTGGTTCGGGCGACCGCGG

Output sequences
Seq1:

TAAGGCGCCAGAACAGCTCACGTCAAACACCCCCGCATAAAGGGAATAAACCAAGTGGACCCTGCTCTAGTGTCTTTG

Seq3:

TTAGACGTCAGAGCAGCTCATGTAAAGCACCCCTAAACAAAGGAAAAAACCAAATGAAATCTGCCCTAATGTCTTTG

Seq2 is absent because it has no match for the primers so it should be discarded.
In seq 3 there is a mismatch in Forward (CGAAAA...) and in Reverse (GTTCGG...).  I don't need the primers to have 100% identity.  85-90% is OK.
Issues:

primer sequences can be anywhere in each sequence, not just in the beginning.

The sequence between the primers is not constant.  it is approximately 80 bases but not exactly 80.

I don't know how to put the primers' identity parameter.

How can I search for R or N nucleotides (IUPAC)? Because I don't know how, I wrote the primers with specific nucleotides and not R or N.

This is the code I have so far, but it has the problems mentioned above:
from Bio import SeqIO

def trim_adaptors(records, adaptorF,adaptorR):

    len_adaptorF = len(adaptorF)
    len_adaptorR = len(adaptorR)

    for record in records:
        indexF = record.seq.find(adaptorF)
        indexR = record.seq.find(adaptorR)
        if indexR == -1  or indexF == -1:
            # adaptor not found, so won't trim
            print(record.id, "no primer/s match")
            continue
        else:
            # trim off the adaptor
            yield record[indexF + len_adaptorF:indexR + 1]

original_reads = SeqIO.parse("trimming_testfas.fas", "fasta")

trimmed_reads = trim_adaptors(original_reads,"CGAGAAGACCCTATGGAGCT" ,"GTTGGGGCGACCGCGG")

count = SeqIO.write(trimmed_reads, "trimmed.fasta", "fasta")

print("Saved %i reads" % count)

This doesn't have to be a script in Python, I am also open to using existing tools.
Thank you a lot in advance, hope you got me :)

The IUPAC ambiguity codes can be thought of as regular expression character classes. R matches any purine, so [AG], Y matches any pyrimidine, so [CT], and N matches anything at all, so . or [ACTG]. So those are easy to handle: instead of  CGAGAAGACCCTRTGRAGCT you can use the regular expression CGAGAAGACCCT[AG]TG[AG]AGCT, and instead of GTTGGGGYGACCNYGG you can use GTTGGGG[CT]GACC.[CT]GG.
Allowing for random mismatches is trickier though and since you're open to using existing tools, I strongly urge you to use something like cutadapt instead of trying to reinvent the wheel. So, install cutadapt on your system (see here for instructions) and then run:
cutadapt -a "CGAGAAGACCCTRTGRAGCT...GTTGGGGYGACCNYGG;max_error_rate=0.15;" --discard-untrimmed -o output.fa file.fa

That assumes your input sequences are in fasta format in the file file.fa and  will create the output file output.fa:
$ cat output.fa
>seq1
TAAGGCGCCAGAACAGCTCACGTCAAACACCCCCGCATAAAGGGAATAAACCAAGTGGACCCTGCTCTAGTGTCTTTG
>seq3
TTAGACGTCAGAGCAGCTCATGTAAAGCACCCCTAAACAAAGGAAAAAACCAAATGAAATCTGCCCTAATGTCTTTG

The options used are -a $primer1..$primer2 which tells cutudapt these primers are linked adapters (https://cutadapt.readthedocs.io/en/stable/guide.html#linked-adapters) and then we also add max_error_rate=0.15; to allow up to 85% identity matches. Cutadapt already knows how to handle IUPAC ambiguity codes so nothing special is needed for those.
I am not sure how this will handle cases where only one of the two sequences is present, you'll need to test for that. I expect that they will be discarded since we're specifying linked primers but I would double check to be sure. In any case, cutadapt is a very powerful tool, so I'm sure there will be a way to handle that in the way you want if you check the documentation.

Thank you a lot! Sequences with only 1 primer (and the 2nd is absent or more mismatches than the given threshold) aren't discarded.  cutadapt cuts the one primer which is present and let the sequence in the dataset which i dont want to.  But generally it works.  Σε ευχαριστώ θερμά! @terdon


text:
Print indel/mismatch in a blastp stream
I am in need of a simple tool that accepts an XML file generated by a BlastP as input and gives me back a file containing information about mismatches/indels/etc. Fortunately such a tool already exists, but it works specifically with XML files generated by BlastN: http://lindenb.github.io/jvarkit/BlastNToSnp.html
I'm trying to adapt the code so that it also reads BlastP files, but after adapting the code it gives me both information about matches and mismatches in the output, and since my knowledge of Java is quite limited I can't identify the problem. Can someone with more expertise in java help me?
package com.github.lindenb.jvarkit.tools.blast;

import java.io.File;
import java.io.PrintWriter;
import java.util.List;

import gov.nih.nlm.ncbi.blast.Hit;
import gov.nih.nlm.ncbi.blast.Hsp;
import gov.nih.nlm.ncbi.blast.Iteration;

import javax.xml.bind.JAXBContext;
import javax.xml.bind.JAXBElement;
import javax.xml.bind.JAXBException;
import javax.xml.bind.Marshaller;
import javax.xml.bind.Unmarshaller;
import javax.xml.namespace.QName;
import javax.xml.stream.XMLEventReader;
import javax.xml.stream.XMLInputFactory;
import javax.xml.stream.XMLResolver;
import javax.xml.stream.XMLStreamException;
import javax.xml.stream.events.XMLEvent;

import htsjdk.samtools.util.CloserUtil;

import com.github.lindenb.jvarkit.io.IOUtils;
import com.github.lindenb.jvarkit.util.bio.AcidNucleics;

import com.beust.jcommander.Parameter;
import com.github.lindenb.jvarkit.util.jcommander.Launcher;
import com.github.lindenb.jvarkit.util.jcommander.Program;
import com.github.lindenb.jvarkit.util.log.Logger;

@Program(name="blastp2snp",
    keywords={"blast","snp"},
    description="print indel/mismatch in a blastp stream",
    biostars={89151,455765},
    modificationDate="20200831",
    creationDate="20131126"
    )
public class BlastPToSnp extends Launcher
{
    private static final Logger LOG = Logger.build(BlastPToSnp.class).make();

    @Parameter(names={"-o","--output"},description=OPT_OUPUT_FILE_OR_STDOUT)
    private File outputFile = null;

    @Parameter(names={"-n","--gapsize"},description="min gap")
    private int minGapSize = 3 ;

    private Unmarshaller unmarshaller;
    private Marshaller marshaller;
    /* force javac to compile */
    @SuppressWarnings("unused")
    private gov.nih.nlm.ncbi.blast.ObjectFactory _ignore_for_javac=null;

    private Iteration peekIteration(final XMLEventReader r) throws XMLStreamException,JAXBException
        {
        while(r.hasNext())
            {
            final XMLEvent evt=r.peek();
            if(!(evt.isStartElement()))
                {
                r.next();
                continue;
                }
            final String name=evt.asStartElement().getName().getLocalPart();
            if(name.equals("BlastOutput_program"))
                {
                r.next();
                String prog=r.getElementText();
                if(!"blastp".equals(prog))
                    {
                    throw new XMLStreamException("Not a blastp input:"+prog);
                    }
                continue;
                }

            if(!name.equals("Iteration"))
                {
                r.next();
                continue;
                }

            return this.unmarshaller.unmarshal(r, Iteration.class).getValue();
            }
        return null;
        }

private void run(
        PrintWriter pw,
        XMLEventReader r
        )
        throws XMLStreamException,JAXBException
    {
     long numIterations=0L;
     long numPerfectMath=0L;
     long numNoHit=0L;

    pw.print("#query");
    pw.print('\t');
    pw.print("hit");
    pw.print('\t');
    pw.print(("hit-index"));
    pw.print('\t');
    pw.print("hsp-index");

    pw.print('\t');
    pw.print("query-POS");
    pw.print('\t');
    pw.print("hit-POS");
    pw.print('\t');
    pw.print("STRAND");

    pw.print('\t');
    pw.print("REF(hit)");
    pw.print('\t');
    pw.print("ALT(query)");
    pw.print('\t');
    pw.print("blast.align_length");
    pw.print('\t');
    pw.print("blast.hit.var");
    pw.print('\t');
    pw.print("blast.query.var");
    pw.print('\t');
    pw.print("blast.mid.var");

    pw.println();
    for(;;)
        {

        Iteration iter1=peekIteration(r);
        if(iter1==null) break;
        ++numIterations;
        if(iter1.getIterationHits().getHit().isEmpty())
            {
            LOG.info("No hit found for "+iter1.getIterationQueryDef());
            ++numNoHit;
            continue;
            }
        boolean found_mismatch=false;
        for(int hit_loop=0; hit_loop <  iter1.getIterationHits().getHit().size();++hit_loop)
            {
            final Hit hit=iter1.getIterationHits().getHit().get(hit_loop);
            for(int  hsp_loop =0;hsp_loop< hit.getHitHsps().getHsp().size();++hsp_loop)
                {
                final Hsp hsp = hit.getHitHsps().getHsp().get(hsp_loop);
                String hspQseq=hsp.getHspQseq();
                String hspMid=hsp.getHspMidline();
                String hspHSeq=hsp.getHspHseq();

                int hsp_query_from=Integer.parseInt(hsp.getHspQueryFrom());
                int hsp_query_to=Integer.parseInt(hsp.getHspQueryTo());
                int hsp_hit_from=Integer.parseInt(hsp.getHspHitFrom());
                int hsp_hit_to=Integer.parseInt(hsp.getHspHitTo());
                final int align_length=Integer.parseInt(hsp.getHspAlignLen());
                final int hit_shift=(hsp_hit_from>hsp_hit_to?-1:1);

                int i=0;
                int query_index = hsp_query_from;
                int hit_index = hsp_hit_from;

                while( i< align_length )
                    {
                    char ch= hspHSeq.charAt(i);
                    char cq= hspQseq.charAt(i);
                    char cm= hspMid.charAt(i);
                    if(cm=='|')
                        {
                        ++i;
                        query_index++;
                        hit_index+=hit_shift;
                        continue;
                        }
                    found_mismatch=true;

                    int j=i+1;
                    for(;;)
                        {
                        int k=hspMid.indexOf(' ', j);
                        if(k==-1 || k-j> minGapSize) break;
                        j=k+1;
                        }

                    String ref=new String(hspHSeq.substring(i, j)).replaceAll("[\\- ]", "");
                    String alt=new String(hspQseq.substring(i, j)).replaceAll("[\\- ]", "");

                    if(hit_shift<0)
                        {
                        StringBuilder sb=new StringBuilder(alt.length());
                        for(int x=alt.length()-1;x>=0;--x)
                            {
                            sb.append(AcidNucleics.complement(alt.charAt(x)));
                            }
                        alt=sb.toString();

                        sb=new StringBuilder(ref.length());
                        for(int x=ref.length()-1;x>=0;--x)
                            {
                            sb.append(AcidNucleics.complement(ref.charAt(x)));
                            }
                        ref=sb.toString();

                        }

                    pw.print(iter1.getIterationQueryDef());
                    pw.print('\t');
                    pw.print(hit.getHitDef());
                    pw.print('\t');
                    pw.print((1+hit_loop));
                    pw.print('\t');
                    pw.print((1+hsp_loop));

                    pw.print('\t');
                    pw.print(query_index);
                    pw.print('\t');
                    pw.print(hit_index);
                    pw.print('\t');
                    pw.print(hit_shift==1?'+':'-');

                    pw.print('\t');
                    pw.print(ref);
                    pw.print('\t');
                    pw.print(alt);

                    pw.print('\t');
                    pw.print(align_length);
                    pw.print('\t');
                    pw.print(hspHSeq.substring(i, j).replace(' ', '.'));
                    pw.print('\t');
                    pw.print(hspQseq.substring(i, j).replace(' ', '.'));
                    pw.print('\t');
                    pw.print(hspMid.substring(i, j).replace(' ', '.'));

                    pw.println();

                    //marshaller.marshal(new JAXBElement<Hsp>(new QName("Hsp"), Hsp.class, hsp), System.out);
                    //pw.println();

                    while(i<j)
                        {
                        ch= hspHSeq.charAt(i);
                        cq= hspQseq.charAt(i);

                        if(ch!='-' && ch!=' ')
                            {
                            hit_index+=hit_shift;
                            }
                        if(cq!='-' && cq!=' ')
                            {
                            query_index++;
                            }
                        ++i;
                        }
                    }
                if(hit_index-hit_shift!=hsp_hit_to)
                    {
                    marshaller.marshal(new JAXBElement<Hsp>(new QName("Hsp"), Hsp.class, hsp),stderr());
                    throw new IllegalStateException("Error expected hit-index= "+hit_index+"-"+hit_shift+" == hsp-hit-to="+hsp_hit_to);
                    }
                if(query_index-1!=hsp_query_to)
                    {
                    marshaller.marshal(new JAXBElement<Hit>(new QName("Hit"), Hit.class, hit),stderr());
                    throw new IllegalStateException("query_index "+query_index+"(1 != hsp_query_to:"+hsp_query_to);
                    }
                if(pw.checkError()) break;
                }

            }//end loop Hit

        if(!found_mismatch)
            {
            numPerfectMath++;
            }

        }//end while read
    LOG.info("ITERATIONS : " +numIterations);
    LOG.info("ONLY_PERFECT_MATCH : " +numPerfectMath);
    LOG.info("NO_HIT : " +numNoHit);
    }

    @Override
    public int doWork(final List<String> args) {
        PrintWriter pw=null;
        XMLEventReader rx=null;
        try
            {

            JAXBContext jc = JAXBContext.newInstance("gov.nih.nlm.ncbi.blast");
            this.unmarshaller=jc.createUnmarshaller();
            this.marshaller=jc.createMarshaller();
            this.marshaller.setProperty(Marshaller.JAXB_FORMATTED_OUTPUT,true);
            this.marshaller.setProperty(Marshaller.JAXB_FRAGMENT,true);
            XMLInputFactory xmlInputFactory=XMLInputFactory.newFactory();
            xmlInputFactory.setProperty(XMLInputFactory.IS_NAMESPACE_AWARE, Boolean.FALSE);
            xmlInputFactory.setProperty(XMLInputFactory.IS_COALESCING, Boolean.TRUE);
            xmlInputFactory.setProperty(XMLInputFactory.IS_REPLACING_ENTITY_REFERENCES, Boolean.TRUE);
            xmlInputFactory.setProperty(XMLInputFactory.IS_SUPPORTING_EXTERNAL_ENTITIES, Boolean.FALSE);
            xmlInputFactory.setXMLResolver(new XMLResolver()
                {
                @Override
                public Object resolveEntity(String arg0, String arg1, String arg2,
                        String arg3) throws XMLStreamException
                    {
                    LOG.info("resolveEntity:" +arg0+"/"+arg1+"/"+arg2);
                    return null;
                    }
                });
            final String inputName= oneFileOrNull(args);
            if(inputName==null)
                {
                LOG.info("Reading from stdin");
                rx=xmlInputFactory.createXMLEventReader(stdin());
                }
            else
                {
                LOG.info("Reading from "+inputName);
                rx=xmlInputFactory.createXMLEventReader(IOUtils.openURIForBufferedReading(inputName));
                }
            pw = super.openFileOrStdoutAsPrintWriter(outputFile);
            run(pw,rx);
            pw.flush();
            return RETURN_OK;
            }
        catch(Exception err)
            {
            LOG.error(err);
            return -1;
            }
        finally
            {
            CloserUtil.close(rx);
            CloserUtil.close(pw);
            }
        }

    /**
     * @param args
     */
    public static void main(final String[] args) {
        new BlastPToSnp().instanceMainWithExit(args);

    }

}

source: Stack Exchange

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

sort: votes


text:
Extracting positions from PAF files in order to extract sequences from a Fastq file with Python
I have used Minimap2 to create a paf file by aligning a Fastq file against itself. Now from this Paf file I can see where the reads overlap, and I want to take these positions, and use them to extract the entire sequences between these positions from the paf file so like this:
If read A starts on position 72 and ends on position 2041 and overlaps here with read B which starts on position 190 and ends on position 2741, I would like to extract these from the Fastq and merge into a contig. But I don't really know how to do that.
query_name,query_length,query_start,query_end,strand,target_name,target_length,target_start,target_end,residue_matches,alignment_block_length,mapping_quality,tp,cm,s1,dv,rl
read=73,2677,81,2674,+,read=86140,3729,190,2779,2082,2627,0,S,4,2,0,6
read=73,2677,89,2674,-,read=78530,52349,15309,17890,1977,2623,0,S,4,1,0,6
read=73,2677,95,2674,-,read=8246,38338,32497,35071,1976,2623,0,S,4,1,0,6
read=73,2677,95,2674,-,read=82661,13059,8472,11069,1960,2636,0,S,4,1,0,6
read=73,2677,100,2669,+,read=79559,11152,5925,8498,1928,2617,0,S,4,1,0,6
read=73,2677,89,2674,+,read=75579,16131,10030,12576,1924,2609,0,S,4,1,0,6
read=73,2677,100,2669,+,read=92760,9659,7065,9617,1904,2613,0,S,4,1,0,6
read=73,2677,89,2674,+,read=95665,81160,4544,7111,1866,2624,0,S,4,1,0,6
read=73,2677,81,2660,+,read=75036,31663,27966,30529,1841,2619,0,S,4,1,0,6

This is a snippet from the PAF file. It's comma delimitted, and below are some of the reads from the Fastq file, the reads vary a lot in length.
@9a65dba6-05c9-49d7-ac0e-2ede30be6036 runid=7f8eec05f1beed68c29c091f78bf84807beb3275 read=39 ch=1736 start_time=2021-01-27T22:53:24Z flow_cell_id=PAG49596 protocol_group_id=210105Chi sample_id=D21-130 barcode=barcode01
AGTATACTTCGTTTCAGTTACGTATTGCTAAGGTTAACACAAAGACACCGACAACTTTCTTCAGCACCTGGTGCTGAAGAAAGTTATCGGTATTTCTTTGTGATCAGCCGTAG
+
%))(&&''&((AA8>33??7++)'*+$$55(1@BFFF>?>=>684.-8:2)&.,11168=>,*/13;;;348:A>;7032)2312%*,&')(001123)),+&,*$%#$%$$$
@edaf3d67-a971-43ed-be5b-3816b985d53f runid=7f8eec05f1beed68c29c091f78bf84807beb3275 read=136 ch=1940 start_time=2021-01-27T22:54:06Z flow_cell_id=PAG49596 protocol_group_id=210105Chi sample_id=D21-130 barcode=barcode01
ATTACTTCGTTCAGTTACGTATTGCTAAGGTTAACACAAAGACACCACGACAACTTTCTTCGCACTATGGTGTTCGCTATTATCAGATTCCTGATGAAGAAGAACTAATGAAATTGGTTTCTATTATTTATTGACTTCGTTCTATAATAGAGCAAAGTTTTCAGCTGGAGAAAATTGGTGAAGAAAAAGCTGAGATTACTGCAAAAGCAACTTCTCTCGTGGTAAATCTGTCACCTCTGCAGAGCATTATCTGAAAGAATGAGGAACTTTCCTGAAGTTATTGTTCAATATGGTTGCTCTTCAAAATTGCTAAGCCAGAACAGCAAGATTGATAACATTCATGCATTAGAAGACGCTCTTTATAGTAAGTTGGTTGGCGTGGTACAGTAGACTGTATTAGTGAATATGATGGAGAATTAGCAGTTATTGACTTCAAGACTTCTAAAGCACCTAAGACCGAAAATGGATTGATGGTTATTTTGCAGCAATGGCATGTATTTATTCATTAGTTCATTGGCATGTGCTAAGAAACTTGTTGTATAATAATGGCATGTGAAGATGGAACCTGTAAAGTTTATGAAGGCGTAATAGAAACATCAACATCAGCAAATTACTCGTACATTACATTAGAAACTTTTTGGATTTTCAACTAGAATTAAATGGAAAATGAATTAACCAAAGCACTAGATAAAAAATTCATGAATGCTGCCAAAAAGTTTTCTCTTGAAGTAGAGAAAATTGTATTGGAAGAAAAAAGTGAATTATATTGATGCTATTGTTCTATTTTGTGGAAAAATAGTATTGAAGTGGACTCTGTTACTAAGTTGATTTCTAAACCTTTAAAGGAAAAATTAAAGGCAGATACTGTAATTTGAACTTCATGAAGAACCTCGTGCCGACTACCCATATAAATACTAGTATCTAGGAGGAATTATGTCTGATTTTTTGAATCTGAACAAGTTCCAAGAATCTATAAAGACATCAACAGATGCAAGAAAGAGATATATTCTAAAGTCTTGTTTGCTTTTTGACATTGCCTCAAAGAAAATGGATCATGTTGATTTATTAGAACAACTTCTAGAGAAGCAAAAGAAACTTCTACATGAGATTGGCTTCTGACGATCCTCGTGCTAAAGATATGATGGAGCAAATTCATCAATCTGCTCAGTTGATTGGATTCCCCTAAGGATTCATCACCTGATGTTTTATTTGAGAATATGAACCTTGGTCTTATTAGATTTGCTTGACATGGGAGCACTGTCGCCCTATAATAGACCAGTACACAAAAGCCCAAATCTAATTACAGAAATTTGGAAATCTATGTCTTTTGCATCACTAAAAAAACAGTCCTCTCTTGGATCTCTAACCGCAAAGTTGGTTAAAGAAGTAGAAGACCAATAAAGCAGGTAACGGGCTGATGAGCGTCTTTGACCAGAGAGTTGATAAAAAGCATGGAGCAGTCTGTTGCAGTTATTCGTTTCTTCCAGCACCAGATAAGAAGGATTCTTCCTTGGGCAAAACTTTATTCTCATGCATTTCAAGGATGAGTGGATGGTACGTGAAAATAACTACTGTTGGTAGTAAGGATCCTGTTTCAGAGTACAATACGGTCTCATGGATTCTGGTTCAGATGCCAATAAGCAGATCGCTCCTTCAATCAAAGGCGTAAGCTTTTCTTACTATAGCAACATCTATGTTGTAGAGGATATACAAATCCGAATGGAAGGAAAAATATGTATTTACAAATTTGGTAAGAAGATCTTGACAGATCATCGCCAATAATACAACCTGAGTTCAGGGATGAGTCACCAATCAATCCATTGATTTCACAA
+
%%&-79888F?9.-90/5400/0(&24105EEH?:39<7++999:8'+23695,,,..77*%&'*.'%).&)882))2*220.7,3CIHA7476935-2734D;:;8+--1;003))%'&$%$$&#%&&&&(-57:;44>:;B=<A4//.%+452,-4('(#$##%&++&,.9;???:CEHJ9GEA@</10.:9<HH2<<;>IE10,45*79:1(3%1,3-*)'''*'''31*++40*$$$$)'689<?A5668:)4<=220#''%$$$$$%%%.++2/%#%$'(#+*)(&#&')#$')*+'>><::<65<:,-:;8)&((--4;9>=ED?<<==;;886)6689<;&&)&%$"$-55;>=+9:29;)'('')#%%%$$&&&%**+112&).739:@=-0.-.290<>>@DABCFGFGABBDEIHD>>02'''(/-6977>64;62,+/33/&(#**+,.7?:--07;:74'')4.''((&%$((*100.)%%$##$##$#$%#%#%$$%*$&#$$$$#$%'&%((4-;<79992..++2589=B?;;<;=?@>@DE=FAE>H@?===ADBA=HK@>=;<+(830-1--($*+%$)%&$'#(&%%$$%.5569CH@=>==A@A=CGHJKKNJFEKGB0'(1/HLJA<<'4:7/?>99ABE??@CJNG==>3?EJFCCROOFEGD>EHJLWTQMMJIMNJ?@CFBBC@>=>9:2/.1>CEF?GAFDC;@5::&8(BBIFJHBCEFF;571*II?34=6;=BFFG:55=;??7.+*,56@KLFA@@@C<,(+/0,>E734867526/;<DB<DEDDG?=F:>D@?CI@EB@CMHEF@<@B@A>80,444=AC7781)%%&&%(%$#$&#&''(6-.35''7'%(()&)&()*>74<8<.)038./0.&@BDICA=>KE@<<CC><7<;8747?=>;66<:;9ANFED4%4762/0//3;<8:DB8:=8+;<>:E?@D@/59<<956/13,=88<B6@?=',:8@><89>@;)()+%.*&''())&./0/,%'(&&''(2263'''%07A<;<<C?:?===B<8:JD6;===?A@CB//2/.(%$14333@<1.)31,5643+.:;89=-/9:+&)-./022149892163/+;56>=;,'''2360-..;ABC@@98;<@>E=<?>AA83C>==C?*.*FGE@;;3-.---AD9;<.;89=??;;?HA/$$'81*.3+(%#$%($$#%%$'&$&$%&%*%*,92'')7;8:9;?;9=???8:@76;>=A>>A<;>8=9>;==A9<8?==;(1+).1AC=ECDFF501&*)&#(5258EAB>@&/285>@?3.;@@<:=<A@><;;3009;5+556<AA?78=7C;>?;67+7:;75>;,/92./,'++.//(()/59<@>>CA8,&&8:894%%(.%/01//.0.,+',/'&,/10)%$$(/,&&$.,***#%%%&&$##$$%%##$'&&')26&1331:;A:CD?A=3666:/.300(,*/#'#'+#1;6444-+(+35<;89<7257,@@>1.++DDA4*&(*'*)*7<<<=873)#%%$#$--(''%(*+7968>6.41-7888:>B:;==FDC:6998;;8:+'&$$'+))'#(%+%')$(+2<?<56:3116708;8%445433.474/-1+--5465'&&%$&++,,088:5:;95;/&*::@99=A::;45;?)$$&+,,&$$#$%00--.0(),''(&''''&&&$$%##$%%&$'224(12;.001+$*)(()'%&%##$$&'%$%$#"$))$0+)'8;79:;93=A?;&%136-,,/;.+*3;<>8?>=:<9<<58?98.#$&%

Another thing that is currently confusing me about my Paf file is the target length and query length which seems to make no sense. It's not the length of the overlap nor the length of the entire segment or read. Can anybody shed some light here on the Paf file and how to extract these positions in Python? Please note that I do not want to use external tools for this, only Python and Python libraries.

The default output from minimap2 is tab-separated, not comma-separated, and doesn't have headers. Could you please show how you produced that CSV file (either code or a command line would be helpful)?

Your example CSV output has a query name of "read=73" for all lines,  but the fastq file doesn't have any sequences with that format for read names. It has channel read sequences, but those are different (and also don't match). Could you please add that read into the fastq example (either by using the read uuids in the CSV file, or renaming the sequence IDs in the example fastq file to match the csv file)?

@gringer the read names do have read=NN as the 3rd field, perhaps that is being picked up by the tool? Or are those the "channel reads" you mentioned?

Each flow cell has a fixed number of sequencing channels; reads from each channel are numbered sequentially, and this is the read=NN metadata you see. This is not unique: it's possible to have the same read number from different channels. The channel read sequence should not be used as the read ID; it's better to use the UUID for that (i.e. follow the standard process for read names and use the value after @).

The CSV file was made by an external script, the headers were added that way, so i could understand what was going on, the internal structure of the PAF is still the same.

Please show the command line or script that did the transformation from the PAF file; the read and reference names have been changed from what I would expect from a direct transformation from a PAF file.

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


text:
Does PGS Catalog limit the number of queries to their FTP server?
I am calculating a polygenic risk score (PGS) from the PGS Catalog using the following command:
nextflow run pgscatalog/pgscalc -profile singularity --input {input.samplesheet} --target_build GRCh38 --min_overlap 0.50 --pgs_id {wildcards.PGS_ID} --outdir output/PGS_Catalog_scores/{wildcards.PGS_ID}/nextflow/'

However, I get the following errror:
ERROR ~ Error executing process > 'PGSCATALOG_PGSCALC:PGSCALC:DOWNLOAD_SCOREFILES ([pgs_id:PGS000000, pgp_id:, trait_efo:])'

Caused by:
  Process `PGSCATALOG_PGSCALC:PGSCALC:DOWNLOAD_SCOREFILES ([pgs_id:PGS000000, pgp_id:, trait_efo:])` terminated with an error exit status (1)

Command executed:

  download_scorefiles -i PGS000000                                    -b GRCh38         -o $PWD -v         -c pgsc_calc/1.3.2

  cat <<-END_VERSIONS > versions.yml
  DOWNLOAD_SCOREFILES:
      pgscatalog_utils: $(echo $(python -c 'import pgscatalog_utils; print(pgscatalog_utils.__version__)'))
  END_VERSIONS

Command exit status:
  1

Command output:
  Download failed: [Errno 104] Connection reset by peer

Command error:
  Download failed: [Errno 104] Connection reset by peer
  Traceback (most recent call last):
    File "/usr/local/lib/python3.10/urllib/request.py", line 1565, in ftp_open
      fw = self.connect_ftp(user, passwd, host, port, dirs, req.timeout)
    File "/usr/local/lib/python3.10/urllib/request.py", line 1585, in connect_ftp
      return ftpwrapper(user, passwd, host, port, dirs, timeout,
    File "/usr/local/lib/python3.10/urllib/request.py", line 2406, in __init__
      self.init()
    File "/usr/local/lib/python3.10/urllib/request.py", line 2415, in init
      self.ftp.connect(self.host, self.port, self.timeout)
    File "/usr/local/lib/python3.10/ftplib.py", line 162, in connect
      self.welcome = self.getresp()
    File "/usr/local/lib/python3.10/ftplib.py", line 244, in getresp
      resp = self.getmultiline()
    File "/usr/local/lib/python3.10/ftplib.py", line 230, in getmultiline
      line = self.getline()
    File "/usr/local/lib/python3.10/ftplib.py", line 212, in getline
      line = self.file.readline(self.maxline + 1)
    File "/usr/local/lib/python3.10/socket.py", line 705, in readinto
      return self._sock.recv_into(b)
  ConnectionResetError: [Errno 104] Connection reset by peer

  The above exception was the direct cause of the following exception:

  Traceback (most recent call last):
    File "/venv/lib/python3.10/site-packages/pgscatalog_utils/download/download_scorefile.py", line 81, in _download_ftp
      with closing(request.urlopen(url)) as r:
    File "/usr/local/lib/python3.10/urllib/request.py", line 216, in urlopen
      return opener.open(url, data, timeout)
    File "/usr/local/lib/python3.10/urllib/request.py", line 519, in open
      response = self._open(req, data)
    File "/usr/local/lib/python3.10/urllib/request.py", line 536, in _open
      result = self._call_chain(self.handle_open, protocol, protocol +
    File "/usr/local/lib/python3.10/urllib/request.py", line 496, in _call_chain
      result = func(*args)
    File "/usr/local/lib/python3.10/urllib/request.py", line 1582, in ftp_open
      raise URLError(exp) from exp
  urllib.error.URLError: <urlopen error [Errno 104] Connection reset by peer>

  During handling of the above exception, another exception occurred:

  Traceback (most recent call last):
    File "/venv/bin/download_scorefiles", line 8, in <module>
      sys.exit(download_scorefile())
    File "/venv/lib/python3.10/site-packages/pgscatalog_utils/download/download_scorefile.py", line 66, in download_scorefile
      _download_ftp(url, path)
    File "/venv/lib/python3.10/site-packages/pgscatalog_utils/download/download_scorefile.py", line 88, in _download_ftp
      if '421' in error.reason and retry < max_retries:
  TypeError: argument of type 'ConnectionResetError' is not iterable
  INFO:    Cleaning up image...

Work dir:
  /work/93/2680d519d19d52cf58aabf469963c1

Tip: view the complete command output by changing to the process work dir and entering the command `cat .command.out`

 -- Check '.nextflow.log' file for details

[Thu Apr  4 16:53:26 2024]
Error in rule pgscalc_from_pgs_catalog:
    jobid: 0
    input: input/samplesheet.csv
    output: output/PGS_Catalog_scores/PGS000000/nextflow/score/aggregated_scores.txt.gz
    conda-env: .snakemake/conda/50fcb51d9c875008db6dcc4478dd238d_
    shell:
        nextflow run pgscatalog/pgscalc -profile singularity --input input/samplesheet.csv --target_build GRCh38 --min_overlap 0.50 --pgs_id PGS000000 --outdir output/PGS_Catalog_scores/PGS000000/nextflow/
        (one of the commands exited with non-zero exit code; note that snakemake uses bash strict mode!)

Shutting down, this might take some time.
Exiting because a job execution failed. Look above for error message
srun: error: node02: task 0: Exited with exit code 1

I am running many PGS calculations in parallel. I think the error has to do with the PGS Catalog limiting the number of queries on their FTP server but I am not sure. However, I do not know what the limit would be. Does anyone know how I can avoid this problem? I could limit the number of queries per hour but I would need to know what the limit is to be optimal.
Note, I have rewritten some paths and the PGS ID for privacy reasons.
Thanks!

Is there anything more informative in the full log? I see "Tip: view the complete command output by changing to the process work dir and entering the command cat .command.out". Did you do that? Anything there?


text:
GATK g CNV code Aborts Before Running Main Code
I have been struggling to get the code I wrote for running GATK g CNV on my Linux machine (Steps Here: https://gatk.broadinstitute.org/hc/en-us/articles/360035531152 ) to work. Specifically it seems to abort a run before running the main code (I even removed my DetermineGermlineContigPloidy step to no avail). Below is my code and then the error:

## Paths To Be Set

echo "My paths have been set"

### Set Bam Directory:
#### Make sure to set a seperate directory for the batch you are doing

BAM_DIR=/pathToThe/Batch_I_Am_Running

### Set the path to the reference genome:
reference=/pathToMy/primary_assembly.fa

#model_dir =/path/to/my/GATK-gCNV

### Set the path to the singularity image:
SINGULARITY_CACHEDIR=/path/to/my/GATK
singularity_image=$SINGULARITY_CACHEDIR/gatk_latest.sif

### Set Number Of Bins
#### Normally 100 to 200*

NUM_BINS=200

### Setting the number of threads ?? Like Cores ??

NUM_THREADS=20

#______________________________________________________________________________________________________________________________________________________________

## Main code

echo "I am about to start the main code"

for c in "${CHROMOSOME_NAMES[@]}"; do

  # loop over the directories containing the BAM files
  for dir in "${BAM_DIR}"/*/; do

    # get the name of the directory without the path
    dir_name="$(basename "${dir}")"

    # get the name of the BAM file without the extension
    bam_name="$(basename "${dir}"/*.bam .bam)"

    echo "Processing file: $bam_name"

    echo "I am starting Step 1: Preprocessing"

    # Step 1: Preprocessing
    singularity exec $singularity_image PreprocessIntervals \
  -R ${reference} \
      -I ${dir}/${bam_name}.bam \
      -L targets.interval_list \
      --bin-length 0 \
      -imr OVERLAPPING_ONLY \
      -O targets.preprocessed.interval_list

    echo "I am starting Step 2: AnnotateIntervals with GC content"

    # Step 2: AnnotateIntervals with GC content
    singularity exec $singularity_image AnnotateIntervals \
  -L chr${c}.interval_list \
      -I ${dir}/${bam_name}.bam \
      -R ${reference} \
  -imr OVERLAPPING_ONLY \
  -O chr${c}.annotated.tsv

    echo "I am starting Step 3: CollectReadCounts"

    # Step 3: Collect read counts for the samples
    singularity exec $singularity_image CollectReadCounts \
  -I ${dir}/${bam_name}.bam \
  -L targets.preprocessed.interval_list \
  --interval-merging-rule OVERLAPPING_ONLY \
  --format TSV \
  -O ${dir_name}.${bam_name}.counts.tsv

    echo "I am starting Step 4: FilteringIntervals based on GC-content & read depth"

    # Step 4: FilterIntervals based on GC-content and read depth
    singularity exec $singularity_image FilterIntervals \
  -L chr${c}.interval_list \
      --annotated-intervals chr${c}.annotated.tsv \
  -imr OVERLAPPING_ONLY \
  "${dir}"/*.counts.tsv \
      --max-panel-norm 1.5 \
      --min-panel-norm 0.3 \
      -O chr${c}.filtered.interval_list

    echo "Not doing Step 5 DetermineGermlineContigPloidy for now"

    # Step 5: DetermineGermlineContigPloidy (Mandatory Step)
   # singularity exec $singularity_image DetermineGermlineContigPloidy \
 # --model-parameters ${model_dir}/model_parameters.txt \
     # --contig-ploidy-calls ${model_dir}/contig_ploidy_calls.tsv \
 # --contig-ploidy-priors chr20XY_contig_ploidy_priors.tsv \
 # --output-prefix ${model_dir}/ploidy \
     # --standardized-output \
     # $(find "${output_dir}" -name "*.calls.tsv" -printf '%p ') \

#    for file in ${input_dir}/*.tsv
#    do

    # Get sample ID from filename
#      sample=$(basename ${file} .tsv)
#      # Run DetermineGermlineContigPloidy command
#      singularity exec ${singularity_image} DetermineGermlineContigPloidy \
#        --input ${file} \
#        --contig-ploidy-priors ${contig_ploidy_priors} \
#        --output-prefix ${output_dir}/${sample}_ploidy \
#        --standardized-output
#    done

    echo "Creating shards via bins"

    # Creating shards via bins
    singularity exec $singularity_image IntervalListTools \
  SplitIntervals \
  -R ${reference} \
      -L chr${c}.interval_list \
  -scatter ${NUM_BINS} \
      -o ${dir}/${c}_shard.interval_list

    echo "Starting the Shard Loop that will run GermlineCNVCaller and PostprocessGermlineCNVCalls"

    # Loop over the shards and run GermlineCNVCaller and PostprocessGermlineCNVCalls
    for shard_file in "${dir}/${c}"*_shard.interval_list; do

      echo "Starting Step 6: Run GermlineCNVCaller (Mandatory Step)"

      # Step 6: Run GermlineCNVCaller (Mandatory Step)
      singularity exec $singularity_image GermlineCNVCaller \
    --run-mode COHORT \
    -L "${shard_file}" \
        -I $(find "${dir}" -name "*.tsv" -printf '%p ') \
        --minimum-interval-median-percentile 5.0 \
        --maximum-zeros-in-sample-percentage 40.0 \
        --annotation-cutoffs "QD=2.0" "SOR=2.0" "FS=10.0" "MQ=50.0" \
        --interval-merging-rule OVERLAPPING_ONLY \
        --output-prefix $(basename "${shard_file}" .interval_list)

      echo "Starting Step 7: Run PostprocessGermlineCNVCalls (Mandatory Step)"

      # Step 7: Run PostprocessGermlineCNVCalls (Mandatory Step)
      singularity exec $singularity_image PostprocessGermlineCNVCalls \
    --model-shard-path "${shard_file}" \
        --cnv-call-path $(basename "${shard_file}" .interval_list).calls.tsv \
        --output-dir postprocess-sm \
        --scatter-count ${NUM_THREADS} \
        --num-scatter-plot-columns 8 \
        --max-zeros-in-sample-percentage 40.0 \
        --annotation-cutoffs "QD=2.0" "SOR=2.0" "FS=10.0" "MQ=50.0" \
        --interval-merging-rule OVERLAPPING_ONLY \
        --minimum-contig-length 46709983 \
        --minimum-interval-median-percentile 5.0

    done

  done

done

    #______________________________________________________________________________________________________________________________________________________________

The ".out" file I received after the run aborted is:

My paths have been set
I am about to start the main code

Why would it abort so soon (I changed the number of Gs used and that did not have any effect)?

Please don't dump 168 lines of code on us and expect us to debug it. You need to make this more focused: give us a minimal script that still produces the error. You haven't even put in the minimal effort to remove commented lines of code here! In this case, the error is trivial, but as a general rule, you need to help us help you.

I will try. I am still getting used to coding. Shouldn't commented lines remain?

No worries, we all start somewhere! Just remember the first rule of debugging: print all the things! By far the most common source of bugs is that the developer assumes their variables have one value, but actually the variables have another value. So the first step, always, should be to print out the variables, and to add simple print statements like echo "I got to here!" at different places in your code so you can figure out exactly where the problem occurs.

That is what I was trying to do with the echo statements.

You are using ${CHROMOSOME_NAMES[@]} but you never define that, so when your code gets here:
for c in "${CHROMOSOME_NAMES[@]}"; do
...
done

It doesn't enter the for loop because "${CHROMOSOME_NAMES[@]}" is an empty array. It's the same thing as doing this in a shell session:
$ for c in "${CHROMOSOME_NAMES[@]}"; do echo "This is chromosome $c"; done
$

Just like your script, that will produce no output since there is nothing for it to loop over.

Thank you for this I will try it out seems that I removed my code CHROMOSOME_NAMES=( $(seq 1 22) X Y)


text:
Where can I find examples of published bioinformatics pipeline packages?
I am a software engineer and I am preparing a presentation to aspiring bioinformatics PhDs on how to use best-practice software engineering when publishing code (such as include documentation, modular design, include tests, ...). In particular my presentation will be focused on publishing "pipeline" code, that is, code that is mainly focused on transforming data to a suitable shape for analysis.
I am trying to find good example of published bioinformatics pipelines that I can point students to, but as I am not a bioinformatician I am struggling to do so, so would like your help. It doesn't matter if the published pipeline is super-niche or not very popular so long as you think it is engineered well so that I can look it up.
Specifically the published code should have: adequate documentation, testing methodology, modular design and be easy to install.

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

@Community I added the list of features that the published package should have to the end of the question.

One quite important point IMO, how are you publishing, e.g. GitHub, Pypy, Docker, Conda? Publishing a data pipeline is not the same as placing monolithic code on GitHub that entirely belongs to the author and has no dependencies subject restricted copyright.

This question is too broad to be useful for the StackExchange format (which works best if there's a single, specific correct answer). I recommend that you ask this at a place where more general discussion is easier, like reddit.

The GATK team maintains a GitHub repository with their pipelines: https://github.com/gatk-workflows. They use the Cromwell Workflow Management System, an open source tool created by the Broad Institute.
You are probably most interested in https://github.com/gatk-workflows/gatk4-data-processing, their workflow for aliging fastq files against a reference genome to produce a bam file which can be used for downstream analyses such as variant calling:

gatk4-data-processing
Purpose :
Workflows for processing high-throughput sequencing data for variant discovery with GATK4 and related tools.
processing-for-variant-discovery-gatk4 :
The processing-for-variant-discovery-gatk4 WDL pipeline implements data
pre-processing according to the GATK Best Practices. The workflow
takes as input an unmapped BAM list file (text file containing paths
to unmapped bam files) to perform preprocessing tasks such as mapping,
marking duplicates, and base recalibration. It produces a single BAM
file and its index suitable for variant discovery analysis using tools
such as Haplotypecaller.

If you are starting with FASTQ files visit the seq-format-conversion repository for workflows to convert FASTQs to unmapped BAMS.
The processing-for-variant-discovery-gatk4 provides quick and general processing for sequence data using the latest releases of
GATK. If users are interested in a more elaborate version of this
workflow with quality control tasks and routinely tested for validity
(useful in production environments) then visit the
gatk4-genome-processing-pipeline repository.
The BAM output from processing-for-variant-discovery-gatk4 can be used to perform a variety of other analysis like somatic short variant
discovery, germline short variant discovery, or germline copy number
variant discovery. Visit the GATK Best Practices documentation to
determine what to do next with the BAM files.

I think the nextflow core (https://nf-co.re) pipelines should satisfy most of your requirements. These pipelines are written in nextflow (https://nextflow.io) which allows scaling from personal to hpc infrastructure rather efficiently.
The core pipelines are usually quite well documented, and have a module design. I have not really looked into the testing methodology, but I think the modules themselves are tested. Running them is easy once nextflow is installed on your system.


text:
CNV Kit ` from . import commands ImportError: cannot import name 'commands' from '__main__'`
I am trying to run some code for my colleague in bash
/path/to/cnvkit.py batch /path/to/my/folder/with/bams/*.bam \
  --normal --targets ${bed_file} \
  --fasta path/to/my/resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta \
  --output-reference /path/to/my/CD_BATCH1_reference.cnn \
  --output-dir /path/to/my/Group_1 --scatter

however, I keep getting this truly peculiar error
Traceback (most recent call last):
  File "/path/to/cnvkit.py", line 4, in <module>
    from . import commands
ImportError: cannot import name 'commands' from '__main__' (/path/to/cnvkit.py)

Attempt 1
I checked if the path was correct and multiple times I found that it was correct. It did not change the error.
Attempt 2
Perhaps I need to chmod +x * everything in the CNV kit directory and that will help. It did not change the error.
Why is this error occurring? How do I solve this issue?

Could you download it from conda and try again? Its here, https://docs.conda.io/en/main/miniconda.html
Then,
conda create --name myenv
conda activate myenv
conda install -c bioconda cnvkit
cnvkit batch /path/to/my/folder/with/bams/*.bam --normal --targets ${bed_file} --fasta path/to/my/resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta --output-reference /path/to/my/CD_BATCH1_reference.cnn --output-dir /path/to/my/Group_1 --scatter

The path to the file is not used, the command is simply cnvkit
This avoids all the privilege stuff and ensures you're running a virtual environment to protect your system Python.

To address the comments. "conda" is slow and doesn't ever load.
There is an alternative to conda called mamba, this is fast and designed to address this exact issue. I have never needed to use it so the following is a rough guide only.
Download and install via
conda create --name myenv
conda activate myenv
conda install -c conda-forge mamba

The documentation is here, https://mamba.readthedocs.io/en/latest/user_guide/mamba.html
I've never used it but the idea is ...
mamba create -n mambaenv <list of packages>
mamba activate mambaenv

The list of packages could be left blank for now, then,
mamba install "Python>=3.10.0" numpy pandas pathlib
mamba install cnvkit # SEE NOTE

Can also do
mamba install -n mambaenv cnvkit

NOTE cnvkit exists strictly in -c bioconda. You do not seem to search and install via channels in mamba
You search for a package and its dependencies via
mamba repoquery search cnvkit
mamba repoquery depends -a cnvkit

It is important to place a virtual environment between the system Python and package Python (don't want a system corruption). More to the question, this would solve all the permission problems without changing your system in anyway (kinda important).

I did try this before and I tried it again now but for whatever reason I get Solving environment: failed with initial frozen solve. Retrying with flexible solve. over and over again. It never works though for the command conda install -c bioconda cnvkit.

Okay @Indira thanks. Addressed above

"I have never needed to use it" you lucky so-and-so! We regular joes need to stick to mamba if we need conda environments solved in one human lifetime.

@RamRS this is funny. I suggest the reason is that "Regular Joses" (pronounced Hosays) where I am, don't use conda. "Regular Joes" where you are, all use conda. Thus your server always times out. :-)


text:
How to add bootstrap values to the phylogenetic tree generated by OrthoFinder?
When we run the OrthoFinder analysis tool on a group of genomes to get the orthologues shared by them one of the output files include a folder named 'Species_Tree' that contains a text file named 'SpeciesTree_rooted.txt' . These are its contents (This is an example I have used for the question - 11 genomes of the Metabacillus genus) :

((Mlacus:0.244776,Mmangrovi:0.284609)1:0.0277466,(((Mdong:0.0284157,Midriensis:0.0250137)1:0.0553715,Mindicus:0.0904605)1:0.148105,(Mfastidiosus:0.207986,(Mniabensis:0.142138,((Mcrassostreae:0.0124653,Mlitoralis:0.0139468)1:0.116616,(Mhalosaccharovorans:0.0408989,Mschmidt:0.0422379)1:0.0872204)1:0.0402912)1:0.0849362)1:0.0750784)1:0.0277466);

I used iTOL to construct the tree using this file but I could only see the branch lengths and no way to add any bootstrap values if there had been any, to the tree as seen in most research papers - 99, 98, 100 etc. Any suggestions?

Bootstrapping is easy just use iqtree
iqtree -s file -m GTR+I+G -T 6 -B 1000

This will parallelise across 6 CPUs and results in 1000 bootstraps, whilst GTR is the mutation model, which models invariant sites under a 4 category gamma distribution. This is a robust model.
The file is aligned fasta format (its all data pipeline so I forget if I've shifted the format) - although this can be changed to any input format
The result will be in the treefile and is a professional quality tree, complete with bootstraps.
The following up analysis is a Bayesian tree and is more complicated.
Note You must resolve/delete gaps using iqtree because of the way the authors are calling gaps - this related to your question here - How to use Gblocks for trimming single-copy gene sequences? . Basically columns with gaps should be deleted.

Question, are you sure this is not a Bayesian analysis? It would make sense if it was a Bayesian analysis that is 100% robust (there are robustness scores there). This would make sense for this tree and this data. As far  I am aware iTOL is used for viewing and manipulating trees rather than creating them.

Thank you very much for the answer @M__ . Yes, iTOL is used for viewing trees and not creating them. The tree itself is created by Fasttree which is part of OrthoFinder, but the output is a text file. I did use this command and I got an error : File not found GTR .  And then I included -st DNA in the command but it still did not resolve the issue.  This is what I get -- WARNING: Number of command-line arguments differs from checkpoint
WARNING: Command-line differs from checkpoint!
Reading alignment file Metabacillus_Tree.fa ... Fasta format detected
...

Hi @KirtiKulanthaivel the best thing to do arrange a chat? If ok you will need your data available and access to a terminal and understand bash commands such as more and/or have access to a text editor such as BBEdit. GTR is central to the model, it will not work without it. Also could you install iqtree via conda? If that is not clear I will explain. If ok please specify availability in GMT, early afternoon is preferred. It is better to do this interactively than via question/answer. I will ask questions and ask you to post a small sample of your data and output errors. Let me know.

The bioinformatics chat room is here: https://chat.stackexchange.com/rooms/58858/bioinformatics. .. can alternatively do this via question/answer - the key thing is GTR must be resolved. Something is not right for this error to be thrown and I don't know what it is. The theory as to why is complicated but getting it working is getting it working.

Thanks @M__ how about 15:00 hours GMT?

@KirtiKulanthaivel my apologies, I am not around today, but tomorrow is fine. 15:00 hours cool.

Sure @M__ . See you!

Hi @M__ ! Just wanted you to know that the IQTREE code worked but I had to increase the bootstrap replicates from 100 to 1000 when the terminal read this message - #replicates must be >= 1000
The tree came out well, thank you very much. By the way, when running iqtree, there is a small table that shows the Gap/Ambiguity, Composition and p-value for the aligned SNP files in question (for 11 genomes in my case). It says 0.00 under Gap/Ambiguity and p-value and 'failed' under composition, although the tree is made. Why is that so?

Hi @KirtiKulanthaivel thanks for the feedback. The gap/ambiguity is good because you removed all "gaps" via trimal. The ambiguity, is the sequence is high quality - no ambiguities. The p-value of 'failed' I would need to look up. It sounds like it's attempting to estimate ambiguities/gaps according to a known distribution (e.g. Poison) but it can't because there aren't any. If it was a p-value <0.01 I suspect that would be a problem, 'failed' seems fine.


text:
regex: samtools command to "refine" PacBio IsoSeq data?
I am working with preprocessed data (IsoSeq PacBio) and I cannot understand one of the steps:
samtools view -b -e [rq] >= 0.9 flnc.bam

What exactly was filtered from the bam file?
To obtain flnc they used the standard IsoSeq pipeline (lima -> isoseq3 refine).

The samtools view command you've provided just uses a filter expression to discard records that do not meet the specified expression criteria. In this case, the expression is looking for alignments where the value of the 'rq' (read quality) tag is greater or equal to 0.9. The isoseq tags documentation says that the read quality of the alignment is really the predicted accuracy of the polished isoform. Note that the command produces BAM output, so be sure to redirect stdout or write to a file for example using:
samtools view -b -e '[rq] >= 0.9' -o out.bam flnc.bam

thank you very much!


text:
Why does the subset of rows that match a condition added by the rows that don't match a condition not add up to the total number of rows?
I need to match the external gene names of a table I have called (which has 1,000,000 entries) with their respective ensemble gene IDs.
Here's a small subset of df:

V1
external_gene_name
V3

M0085
A1BG
1

M0085
SNX8
1

and here's a small subset of ensembl_ids:

ensembl_gene_id
external_gene_name

ENSG00000288263
GRK1

ENSG00000291428
BRSK2

I also included a line to see the number of entries in external_gene_name which do not have corresponding ensemble names. Here is my code:
merged_df <- merge(df, ensembl_ids, by = "external_gene_name")
elements_not_in <- df$external_gene_name[!(df$external_gene_name %in% ensembl_ids$external_gene_name)]

But what's confusing me is that when I add the number of rows in merged_df with that of elements_not_in, the total exceeds 1,000,000. Shouldn't the total be 1,000,000? I'm wondering if there's a flaw in my code.

From a brief glance, the gap is not in your code but in your expectations. You're operating based on an assumption of 1:1 mapping between ENSG IDs and HGNC symbols (your external_gene_name attribute). In reality, multiple ENSG IDs will map to the same HGNC symbol so the merge will produce multiple rows per external_gene_name in the merge product. This is why you're seeing a total greater than the sum of number of unique entries when you compare the merge product + the anti-join result (i.e., the elements_not_in dataset) to the initial individual datasets.


text:
Merge fragment reads
I'm writing a function (python) that merges overlapping paired reads of a fragment, and a corresponding reference sequence. There's quite a lot of edge cases (e.g. when the reads disagree on indels) so I wonder if there's some standard library that does that?
Context - the input is pysam AlignedSegment DNA reads. The output is fragment, reference - both strings corresponding to the merged fragment and its corresponding reference.
e.g.
read1:           ACGTACGTACGT
read2:                CGTACGTACGT
expected output: ACGTACGTACGTACGT

read1:           ACGTANGTACGT
read2:               ACGTACGTACGT
expected output: ACGTACGTACGTACGT

this is just an illustration. The idea is to merge the overlapping DNA read sequences into the fragment and its corresponding reference using the information in the AlignedSegment (i.e. query_sequence, reference_start, query_qualities, cigartuples, and get_reference_sequence())
Thanks.

Thanks. But why is this not an assembler? This is pretty much exactly what an assembler does, it assembles reads into contigs. Or, if you already have the reference genome, what a read aligner does.

well, I mean it's not an aligner in the sense that's not aligning reads to a reference but to each other.

OK, so it's an assembler. But then what do you mean by the word "reference" in your question? Are you referring to the consensus sequence you can build by assembling the reads, the inferred genomic source sequence? Also, do you absolutely have to write this from scratch and in a scripting language? Can't you use one of the many existing genome assemblers?

I believe the OP wants to work at the "read-pair level" to obtain contigs of mate pairs and not all of the available reads. If that is the case, I would check make.contigs from the mother suite.

I would use existing high-performance software for this, for example BBmerge from BBmap suite, or FLASH, probably many others exist as well. I see no point in re-writing new software when it already exists, unless this is a learning exercise, or you think you can do better algorithm-wise, but then, with respect, you would probably not be asking here.

I'd suggest looking at similar tools that already exist, e.g. here is one that I hadn't heard of that lives on github. Just glancing at the code, it seems to use a scoring heuristic to decide how to derive a merge point, so it's using that to get around the absence of unambiguous points, probably at some rate of error.
BioPython has some functions that will create consensus sequences for you as well based on a pairwise alignment. (Bio.Align.AlignInfo.SummaryInfo.dumb_consensus, linked in that post.)
As a command line tool pandaseq does the same thing, I'd suggest looking at how they handle edge cases.
There are similar utilities in QIIME, NGMerge, BBMerge, [illumina-utils]- this is a pretty well-trodden path, I think all of those will be illuminating to look at in terms of reusing or writing original code.
I'm guessing it will be easier to merge FASTQs and then align, rather than align and then merge. (Note that these tools don't seem to operate on alignments, because that is trickier).
Alternately, if your real goal is to create a consensus genome sequence from your reference based on your alignments, there are a large number of tools that do that, notably samtools (see e.g. https://www.biostars.org/p/419927/). No new code needed, you just input your BAM file and reference sequence.


text:
how to use BLAST+ in an optimal way
I have a 32-bp query that I want to perform BLAST+. my query represents different RNA viruses and I do expect that my hits get more than 500,000.
but in practice, when I perform it I get around only 700 hits. would you please kindly help me how to increase the efficacy of my BLAST+?
I leave my syntax here.
blastn -task megablast -db nr -outfmt 7 -remote -max_target_seqs 99999999 -perc_identity 80 -query D:/Query.fas  -out D:/result.txt

Can you explain why you expect many queries? More importantly, what is the biological question you are trying to answer here? Are you looking for distant homologs? Close homologs? Why did you choose to use megablast? That's an odd choice for such a short sequence.

Thanks for your comment.
I expect many hits because there are many hits for my query as I know my query is common among RNA viruses, it is such a conservative seq.
Yes, I want to extract all sequences in the NR database that are similar to my query.
I think megablst was a bad choice, because when I choose -task blastn my hits were significantly improved (+500,000), but still I encounter an error by blast+ stating "search cannot be formatted" resulting in finding just some of those hits.

Please [edit] your question and explain this. Megablast is only useful for very similar sequences. And long ones, usually. It doesn't really make sense to use it for a 32bp query, that will be smaller than the smallest thing it can find.  So please [edit] and explain what you are trying to do. You don't want similar sequences, presumably you are looking for homologs, right? I mean explain the biological question you are trying to answer, so we know how to help.

Megablast is designed to find very similar sequences. It is highly specific but has a very low sensitivity. This is because it has a very large word size. By default, I think it will use a word size of 28 nt. That's almost your entire sequence, so if your target sequence doesn't have 28 contiguous nucleotides that are identical to 28 contiguous nucleotides of your query, no hit will be found.
If you are looking for distant homologous sequences, you should use simple blastn which defaults to a word size of 11. That will still likely be too small, you might want to change it to something lower (your blast tool should have an option for that), but it will give you more results than megablast.
Finally, if these 32nt are translated into protein, you can consider using tblastn instead.


text:
How can I limit blast results to exclude a specific taxonomic ID?
I would like to perform blast+, using the command line, here is the code:
-task blastn -db nr -outfmt 7 -query myQuery.fas -out result.txt -remote.

I have to exclude SARS-CoV-2 hits from my results, so I must include -negative_taxids mentioning SARS-CoV-2 NCBI taxid number. However, the problem is, -negative_taxids is NOT compatible with -remote command and my command line encounters an error.
I have to use -remote as I cannot download entire +400 GB NR database to run it locally.
Could you please kindly help me out with this contradiction in my work?
My output txt file for this command line is:
# Fields: query acc.ver, subject acc.ver, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
# 600483 hits found
AstV-VA1    MN464146.1  100.000 32  0   0   1   32  6362    6393    7.14e-06    59.0
AstV-VA1    MN148428.1  100.000 32  0   0   1   32  6685    6716    7.14e-06    59.0
AstV-VA1    MK671313.1  100.000 32  0   0   1   32  6682    6713    7.14e-06    59.0
AstV-VA1    MK671312.1  100.000 32  0   0   1   32  6682    6713    7.14e-06    59.0
AstV-VA1    MK671311.1  100.000 32  0   0   1   32  6682    6713    7.14e-06    59.0

Have you considered running your query against all of nr and then removing the SARS hits from the results instead?

Thanks terdon for your constant and positive contribution to my questions.
Actually, I have already done it and I got +600,000 hits but the point is the txt result file contains only the following information:
"query", "acc.ver", "% identity", "alignment length", "mismatches", "gap opens", "q. start", "q. end", "s. start", "s. end", "Evalue", "bit score"
so, there are no organism and taxid columns in the result. One possible solution can be using R to download the annotation to add the results automatically, but it has two disadvantages:

downloading annotation of +600,000 hits will take around 200 hours (more than a week)
not all hits automatically yield annotations.
Therefore, I think filtering the results from the beginning can be a better option.

No, no, there is no reason to download entire annotations! All you need to do is extract the IDs of the target sequences and then get the taxids for those sequences. If you [edit] your question and show us a few examples of your output, we should be able to help, but it will just be a simple entrez query whose input is the list of sequence identifiers and output will be the associated taxids.

What do you want to do with these BLAST hits?

This isn't particularly elegant, but it should work. First, extract the accessions into a file. Next, get the taxid for each of the accessions. Once you have that, you can identify the accessions that are for SARS-CoV-2 sequences, and remove them from your output. Something like this:

Extract the accessions into a new file
awk '!/^#/{print $2}' blast.out > names

Map those accessions to their taxids.
while read -r accession; do
   printf '%s\t%s\n' "$accession"  \
       $(curl -s "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=nuccore&id=${accession}&rettype=summary&retmode=text" |
          awk -F'[<>]' '/Name="TaxId"/{print $3}');
done < names > taxids

With your example input, that produces:
$ while read -r accession; do
>    printf '%s\t%s\n' "$accession"  \
>        $(curl -s "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=nuccore&id=${accession}&rettype=summary&retmode=text" |
>           awk -F'[<>]' '/Name="TaxId"/{print $3}');
> done < names
MN464146.1   1435485
MN148428.1   2661902
MK671313.1   1239566
MK671312.1   1239566
MK671311.1   1239566

Now, this will probably be way too slow, doing them one at a time. So if you have many IDs, try and combine multiple queries into one. Like this:
perl -ne 'chomp; print $. % 50 == 0 ? "$_\n" : "$_," ' names |
 | sed 's/,$/\n/' |
    while read accession; do
      curl -s "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=nuccore&id=${accession}&rettype=summary&retmode=text" |
        awk -F'[<>]' '/Name="Caption"/{ acc=$3 }
                   /Name="TaxId"/{printf "%s\t%s\n", acc, $3}'
    done  < names > taxids

Here, the first Perl command will read through the list of accessions and print the accession followed by a comma and, once every 50 lines ($. % 50 == 0), it will instead print it with a newline. The sed command will replace a trailing comma with a newline for the very last line of output where, unless the number of accessions is an exact multiple of 50, you will have a line ending with a comma instead of a newline.
The result is a list of 50 comma-separated identifiers per line (fewer, on the last line). This is then fed to the curl command and then a naive awk that parses it and prints out the accession and the taxid.
Important: this approach loses the version information. So MK671311.1 becomes MK671311. This shouldn't be a problem, and we will take it into account later, but it should be understood. This is what the output looks like:
MN464146    1435485
MN148428    2661902
MK671313    1239566
MK671312    1239566
MK671311    1239566

Remove the SARS-CoV-2 sequences.
The taxid for SARS-CoV-2 is 2697049, so you can get the list of its accessions with:
awk -F'\t' '$2=="2697049"{print $1}' taxids

You can now filter your file to remove them using:
awk '(NR==FNR){
       if($2=="2697049"){
         sars[$1]
       }
       next
      }
      {
        acc=$2;
        sub(/\.[0-9]*$/, "", acc);
        if(acc in sars){ next }
          print
      }' taxids blast.out > blast.filtered.out

This awk script will read the first input file (NR is the current input line number while FNR the line number of the file being read; the two are only equal while reading the first file) and if its second field is the desired taxid, it saves the first field, the accession, in the array sars. The next ensures we don't do anything else while reading the first file.
When we move to the second file, we make a copy of the second field, the accession, and save it as acc, and then we remove the version information ( sub(/\.[0-9]*$/, "", acc). If the resulting accession is found in the sars array, we skip the line with next and if not, we print.
This assumes you are using the faster, batch retrieval I gave above, the one that removes the versioning information. If you instead use the very slow first one, you can filter with
 awk '(NR==FNR){ if($2=="2697049"){ sars[$1]} next }
      { if($2 in sars){next} print}' taxids blast.out > blast.filtered.out

That should give you a filtered blast output file that no loner has sequences from taxid 2697049. However, please make sure to test this, count the number of input and output lines and make sure everything makes sense. I don't have access to your actual data and while this does work with the example you provided, I cannot guarantee it will work out of the box for your real data. As always, you need to be careful!

@FarzadBeikpour kindly consider upvoting. This is a good answer and a lot of work. If you can't get -negative_taxids working I'd use this approach

I would personally look at is the version number of the blast you are using.
Bugs around blast+ versions >= 2.9.0 (April 2019) for -taxids and I would assume -negative_taxids flag should have been fixed as described in the the Blast+ release notes here.

Fix taxID filtering combined with mask-based alias BLAST databases.

Bug fix descriptions are usually short, hence the lack of further detail.
Blast taxids can be implemented via lists here https://www.ncbi.nlm.nih.gov/books/NBK569846/


text:
error in downloading the annotation of sequences on NCBI
I have developed an R function to download the "organism" name of the numerous NCBI accession numbers which are already saved in an object.
The function:
get_organism_names <- function(accession_numbers) {
  organism_names <- sapply(accession_numbers, function(accession) {
    record <- entrez_summary(db = "nucleotide", id = accession)
    return(record$organism)
  }, USE.NAMES = FALSE)
  names(organism_names) <- accession_numbers
  return(organism_names)
}

Although it works for some accession numbers, for others it encounters an error stating for example:
Warning: Esummary includes error message: Invalid uid KU253616.1 at position= 0Error: No esummary records found in file

How can I either solve this error or ask R to skip the problematic accession and keep downloading the rest?

Have you tried running it in a loop with a single accession each time instead of giving it all of them?

Also, please give us all of the code. Presumably, you are loading some libraries, which ones? We need example code that we can use to reproduce the problem. See https://stackoverflow.com/help/minimal-reproducible-example. We need example input, and then code that we can run locally that reproduces the error.

You could use the tryCatch function to handle inputs that throw an error. A good example on how to write them can be found here. Though the message indicates it's a warning, so does the function stop running when you get this message?


text:
How can I convert a BED file to GTF/GFF with gene_ids?
Given a .bed (BED12), how can I convert it to GTF/GFF formats with gene_id attributes? What is the fastest way or available tools to do it?
For example, given an input like this:
chr27 17266469 17281218 ENST00000541931.8 1000 + 17266469 17281218 0,0,200 2 103,74, 0,14675,

How can I convert it to:
chr27 bed2gtf gene 17266470 17285418 . + . gene_id "ENSG00000151743";
chr27 bed2gtf transcript 17266470 17281218 . + . gene_id "ENSG00000151743"; transcript_id "ENST00000541931.8";
chr27 bed2gtf exon 17266470 17266572 . + . gene_id "ENSG00000151743"; transcript_id "ENST00000541931.8"; exon_number "1"; exon_id "ENST00000541931.8.1";

Here I provide an ordered list of options (note that I am the author of bed2gtf and bed2gff):
bed2gtf
A high-performance BED-to-GTF converter written in Rust from https://github.com/alejandrogzi/bed2gtf.
Usage: bed2gtf[EXE] --bed/-b <BED> --isoforms/-i <ISOFORMS> --output/-o <OUTPUT> --threads/-t <THREADS>

where:
--bed <BED>: a .bed file
--isoforms <ISOFORMS>: a tab-delimited file
--output <OUTPUT>: path to output file (*.gtf)

The isoforms file specification:

a tab-delimited .txt/.tsv/.csv/... file with genes/isoforms (all the transcripts in .bed file should appear in the isoforms file):

> cat isoforms.txt

ENSG00000198888 ENST00000361390
ENSG00000198763 ENST00000361453
ENSG00000198804 ENST00000361624
ENSG00000188868 ENST00000595977

Converts

Homo sapiens GRCh38 GENCODE 44 (252,835 transcripts) in 3.25 seconds.
Mus musculus GRCm39 GENCODE 44 (149,547 transcritps) in 1.99 seconds.
Canis lupus familiaris ROS_Cfam_1.0 Ensembl 110 (55,335 transcripts) in 1.20 seconds.
Gallus galus bGalGal1 Ensembl 110 (72,689 transcripts) in 1.36 seconds.

bed2gff
A Rust BED-to-GFF3 translator that runs in parallel from https://github.com/alejandrogzi/bed2gff.
Usage: bed2gff[EXE] --bed/-b <BED> --isoforms/-i <ISOFORMS> --output/-o <OUTPUT> --threads/-t <THREADS>

where:
--bed <BED>: a .bed file
--isoforms <ISOFORMS>: a tab-delimited file
--output <OUTPUT>: path to output file (*.gff)

The isoforms file specification:

a tab-delimited .txt/.tsv/.csv/... file with genes/isoforms (all the transcripts in .bed file should appear in the isoforms file):

> cat isoforms.txt

ENSG00000198888 ENST00000361390
ENSG00000198763 ENST00000361453
ENSG00000198804 ENST00000361624
ENSG00000188868 ENST00000595977

Convert

Homo sapiens GRCh38 GENCODE 44 (252,835 transcripts) in 4.16 seconds.
Mus musculus GRCm39 GENCODE 44 (149,547 transcritps) in 2.15 seconds.
Canis lupus familiaris ROS_Cfam_1.0 Ensembl 110 (55,335 transcripts) in 1.30 seconds.
Gallus gallus bGalGal1 Ensembl 110 (72,689 transcripts) in 1.51 seconds.

bedToGenePred + genePredToGtf + refTable
UCSC offers a fast way to convert BED into GTF files through KentUtils or specific binaries using:
bedToGenePred in.bed /dev/stdout | genePredToGtf file /dev/stdin out.gtf

You can install these tools with bioconda, or download them here. The gene_id is only achieved when using refTables (a format specified in UCSC's web browser), you can see a more elaborate answer here Obtaining Ucsc Tables Via Ftp And Converting Them To Proper Gff3 Via Genepredtogtf?.
Other options
Other scripts/tools That DO NOT produce a complete GTF file (lacking gene_id attributes) are:

gtf2bed

gtf2bed < foo.gtf | sort-bed - > foo.bed
awk '{print $1"\t"$7"\t"$8"\t"($2+1)"\t"$3"\t"$5"\t"$6"\t"$9"\t"(substr($0, index($0,$10)))}' foo.bed > foo_from_gtf2bed.gtf

-kscript
from https://github.com/holgerbrandl/kscript:
kscript https://git.io/vbJ4B my.bed > my.gtf

pfurio/bed2gtf

from https://github.com/pfurio/bed2gtf:
python bed2gtf [options] <mandatory>

AGAT

AGAT
Considering only the options that produce gene_ids attributes, bed2gtf and bed2gff are faster by ~3-4 seconds than UCSC's C binaries. More detailed instructions of this tools are explained in the sources linked.

Thank you. These are great resources.


text:
Find protein in DNA sequences with arbitrary encoding and possible frameshifts
I have a dataset with lots of DNA sequences (~10 M seqs, each ~6 KBp long).
I want to detect presence of a fixed given protein (~600 Bp / 200 aa) within each sequence and, if detected, obtain the frame and position of the protein within the sequence. However, the protein may be coded arbitrarily and differently in each sequence. In addition, I need to account for frameshifts induced by sequencing errors.
I've tried to use tfasty36 from the FASTA36 package. This is basically a tblastn equivalent (protein query search against DNA database in all frames). Due to the statistical filtering built into these tools, I've come to the conclusion that tfasty36 and tblastn are not the right tools for this job.
tfasty36 didn't work because it requires fiddling with statistical tolerances to show all matches. Also, the matches are quite fragmented (I guess, it uses heuristics instead of dynamic programming). For tblastn, building the database seems to be unnecessary (given that it may have the same problems as tfasty36 afterwards).
I've looked into alignment tools, but couldn't find any that can align a protein sequence to a DNA sequence while taking frameshifts into account.
How can I solve this problem?

If you need to model introns and splicing, then you should use a tool that is designed for that. One of the limitations of blast and its cousins is that they only look at sequence similarity and do not attempt to join HSPs into possible transcript models.
I haven't worked on this in a long time, but back when I did, the programs I used were:

exonerate

Exonerate is a generic tool for pairwise sequence comparison. It allows you to align sequences using a many alignment models, either exhaustive dynamic programming or a variety of heuristics.

In this case, you would use your protein sequence as the query and the set of  DNA sequences as the target and try to align using the protein2genome model. Something like this:
  exonerate --model protein2genome protein.pep dna.fa > out

This will result in (many) alignments of your query protein agains the various genomic sequences. Here's an example I had handy:
$ exonerate --model protein2genome -n 1 frog.histone.pep c.elegans.chrI.fa
Command line: [exonerate --model protein2genome -n 1 frog.histone.pep c.elegans.chrI.fa]
Hostname: [oregano]

C4 Alignment:
------------
         Query: frog
        Target: CM041006.1 Caenorhabditis elegans strain DL226 chromosome I, whole genome shotgun sequence:[revcomp]
         Model: protein2genome:local
     Raw score: 732
   Query range: 32 -> 300
  Target range: 3039620 -> 3034073

      33 : LeuGlnLysLysLeuGluGluLeuGluLeuAspGluGlnGlnArgLysArgLeuGl :      51
           |||:!!:!!:!!|||:!!:!!|||! !|||   |||! !|||! !:!!||||||..
           LeuGluGluGlnLeuLysLysLeuGlyLeuThrGluProGlnThrGlnArgLeuSe
 3039620 : CTAGAAGAACAGCTGAAAAAGCTCGGGCTCACGGAGCCACAAACTCAACGGCTCAG : 3039566

      52 : uAlaPheLeuThrGlnLysGlnLysValGlyGluLeuLysAspAspAspPheGluL :      70
           .! !||||||  !   |||:!!   :!!  !||||||!..!!:|||   .!.:!!!
           rGluPheLeuGlnValLysGluGlyIleLysGluLeuSerGluAspMetLeuGlnT
 3039565 : CGAATTCTTGCAGGTCAAGGAGGGCATCAAGGAGTTGAGTGAGGATATGCTACAGA : 3039509

      71 : ysValSerGluLeuGlyAlaGlyAsnGlyGlyValValPheLys  >>>> Targe :      85
             ! !.!!|||||||||  !||||||||||||||||||  !|||
           hrGluGlyGluLeuGlyHisGlyAsnGlyGlyValValAsnLys++
 3039508 : CTGAAGGAGAACTTGGCCATGGAAATGGAGGAGTTGTTAATAAGgt.......... : 3039462

      86 : t Intron 1 >>>>  ValSerHisLysProThrSerLeuIleMetAlaArgLys :      97
           55 bp              !  !|||::!  !|||.!!:!!|||||||||||||||
                          ++CysValHisArgLysThrGlyValIleMetAlaArgLys
 3039461 : ...............agTGTGTTCATCGAAAAACCGGTGTAATTATGGCTCGGAAG : 3039373

      98 : LeuIleHisLeuGluIleLysProAlaIleArgAsnGlnIleIleArgGluLeuGl :     116
           |||:!!||||||||||||||||||:!!:!!|||.!.||||||:!!!:!||||||
           LeuValHisLeuGluIleLysProSerValArgGlnGlnIleValLysGluLeuAl
 3039372 : CTTGTACACCTCGAGATCAAGCCGTCAGTCCGACAACAGATCGTCAAAGAGCTCGC : 3039316

     117 : nValLeuHisGluCysAsnSerProTyrIleValGlyPheTyrGlyAlaPheTyrS :     135
            |||||||||:!!||||||||||||!:!||||||||||||||||||||||||   .
           aValLeuHisLysCysAsnSerProPheIleValGlyPheTyrGlyAlaPheValA
 3039315 : CGTACTTCACAAGTGCAATAGTCCATTCATTGTTGGATTCTACGGGGCTTTTGTGG : 3039259

     136 : erAspGlyGluIleSerIleCysMetGluHisMetAspGlyGlySerLeuAspGln :     153
           .!:!!..!!!:||||||||||||||||||:!!|||||||||  !|||||||||  !
           spAsnAsnAspIleSerIleCysMetGluTyrMetAspGlyLeuSerLeuAspIle
 3039258 : ACAATAACGATATATCAATTTGTATGGAGTATATGGATGGGCTTTCACTTGATATA : 3039205

     154 : ValLeuLysLys  >>>> Target Intron 2 >>>>  AlaGlyLysIlePro :     162
           ||||||||||||           3837 bp           !.!|||::::!!|||
           ValLeuLysLys++                         ++ValGlyArgLeuPro
 3039204 : GTTTTGAAAAAGgt.........................agGTTGGCCGCCTCCCG : 3035341

     163 : GluLysIleLeuGlyLysValSerIleAlaValIleLysGlyLeuThrTyrLeuAr :     181
           ||||||.!!:!!|||::!:!!|||:!!||||||:!!:::|||||||||||||||!:
           GluLysPheValGlyArgIleSerValAlaValValArgGlyLeuThrTyrLeuLy
 3035340 : GAAAAGTTTGTCGGACGAATCTCGGTTGCAGTTGTCCGTGGATTGACGTACTTGAA : 3035284

     182 : gGluLysHisLysIleMetHisArg{A}  >>>> Target Intron 3 >>>>  :     190
           !!!::!!  !||||||:!!||||||{|}            50 bp
           sAspGluIleLysIleLeuHisArg{A}++                         +
 3035283 : AGATGAAATTAAGATTTTGCATAGA{G}gt.........................a : 3035208

     191 :  {sp}ValLysProSerAsnIleLeuValAsnSerArgGlyGluIle--LysLeuC :     207
            {||}|||||||||||||||!!:||||||||||||!.!|||||||||##   ||||
           +{sp}ValLysProSerAsnMetLeuValAsnSerAsnGlyGluIle##---LeuC
 3035207 : g{AT}GTTAAACCATCGAATATGCTTGTCAATAGTAACGGAGAAATTAA---TTAT : 3035157

     208 : ysAspPheGlyValSerGlyGlnLeuIleAspSerMetAlaAsnSerPheValGly :     225
           ||||||||||||||||||||..!|||||||||||||||||||||||||||||||||
           ysAspPheGlyValSerGlyMetLeuIleAspSerMetAlaAsnSerPheValGly
 3035156 : GCGATTTTGGAGTCTCTGGAATGTTGATTGATTCGATGGCCAACTCATTTGTCGGA : 3035103

     226 : ThrArgSerTyrMetSer  >>>> Target Intron 4 >>>>  ProGluArg :     234
           |||||||||||||||:!!            799 bp           |||||||||
           ThrArgSerTyrMetAla++                         ++ProGluArg
 3035102 : ACGAGGAGTTATATGGCGgt.........................agCCCGAACGA : 3034277

     235 : LeuGlnGlyThrHisTyrSerValGlnSerAspIleTrpSerMetGlyLeuSerLe :     253
           |||  !|||:!!||||||!!!:!!..!|||||||||||||||.!.|||||||||||
           LeuThrGlySerHisTyrThrIleSerSerAspIleTrpSerPheGlyLeuSerLe
 3034276 : CTCACAGGATCCCACTATACAATTTCATCGGATATCTGGTCATTTGGATTATCACT : 3034220

     254 : uValGluMetAlaIleGlyArgTyrProIleProProProAspAlaLysGluLeuG :     272
           |||||||:!:  !|||||||||||||||:!:||| !!|||..!  !  !|||! !!
           uValGluLeuLeuIleGlyArgTyrProValProAlaProSerGlnAlaGluTyrA
 3034219 : TGTCGAACTTCTTATTGGACGATATCCAGTGCCTGCACCGTCTCAAGCAGAATATG : 3034163

     273 : luLeuIlePheGlyCysSerVal<->GluGlyAspProAlaSerSerGluLeuAla :     289
               !!!:|||..!   :!!! !   |||  !!!:! !|||..!|||     !! !
           laThrMetPheAsnValAlaGluAsnGluIleGluLeuAlaAspSerLeuGluGlu
 3034162 : CTACAATGTTTAATGTGGCTGAAAATGAGATTGAACTGGCTGATAGTCTTGAAGAG : 3034109

     290 : ProArgProArgProProGlyArgProIleSer :     300
           |||!.!  !!.!||||||.!!!.!|||   |||
           ProAsnTyrHisProProSerAsnProAlaSer
 3034108 : CCCAACTATCATCCGCCGAGCAATCCGGCGTCA : 3034074

vulgar: frog 32 300 . CM041006.1 3039620 3034073 - 732 M 52 156 5 0 2 I 0 51 3 0 2 M 73 219 5 0 2 I 0 3833 3 0 2 M 32 96 S 0 1 5 0 2 I 0 46 3 0 2 S 1 2 M 14 42 F 0 2 G 1 0 M 26 78 5 0 2 I 0 795 3 0 2 M 48 144 G 0 3 M 21 63
-- completed exonerate analysis

In addition, exonerate lets you control the output you want. Check out the --ryo (roll your own output) option in man exonerate.

genewise

GeneWise compares a protein sequence to a genomic DNA sequence, allowing for introns and frameshifting errors.

Unlike exonerate, genewise only does protein to DNA alignment, but like exonerate it models introns and splicing. 20 years ago, when I did this every day, my impression was that genewise was slower but slightly more sensitive than the default exonerate options (exonerate can be tweaked to be made more sensitive). So I would usually try exonerate first and fall back to genewise if it failed to find a good match.


text:
Downloading a list of all RCSB-PDB entries
I'm trying to download a list of all of the non-redundant protein chain IDs in the RSCB-PDB (not the actual FASTA entries, just the IDs), filtered by some % identity.
Does anyone know how to do this?  I haven't been able to figure it out with the advanced search function .

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.

% identity to what? To each other? To a specific query sequence?

Please add some additional context to your question. Why dod you want to download this list?

@terdon I believe he needs the ones not identical AKA different ones ' non-redundant' as he says

So the ones that do not have more than X% identity to any others in the database, @pippo1980? OK, so what identity? And won't all of them be non-redundant? I mean, if they are identical, they wouldn't have a different ID, would they? (I know, my deep, deep ignorance of the protein structure field is probably evident here)

@terdon I dont know a lot about it either but I guess you are right the OP should specify a number and whath the number stands for it coukd be both sequence identity or structure identity but probably for the latter we are talking about mainchain Root Mean Square Deviation (RMSD),

but maybe is just this : https://www.rcsb.org/news/feature/61b8e5e869569d045d35a9a3 Remove Redundancy in Search Results

@terdon I'm a computational biologist with a background in biology, not computers. I did my PhD work on gene prediction and comparative evolution, my post-doc on protein-protein interaction networks and I have been working in NGS analysis, specifically human variant calling and annotation, for the past several years. Jack of many trades and master of none. .. I know, my deep, deep ignorance of the protein structure field is probably evident here ... YES SURE .... ;-)

Did you try RCSB sequence similarity search?
These two announcement highlight the options you ask about:

Improved Searching with New Sequence Identity Cutoff (2017)
Remove Redundancy in Search Results (2015)

Alternatively, if you'd like to calculate identity on your own, you can download sequences of all the PDB entries from:
https://files.wwpdb.org/pub/pdb/derived_data/pdb_seqres.txt
(270MB).

Thank you!  The https://files.wwpdb.org/pub/pdb/derived_data/pdb_seqres.txt link was what I needed.  Any way to filter them by identity from the text file?

I suppose you'd need to do pairwise sequence alignment and calculate identity.


text:
Sequence Alignment for sequences with the same length
I am doing research on a new method of optimizing sequence alignment process (Needleman - Wunsch algorithm) but the idea would only work with sequences that have the same length. I am wondering if that will have any real-life utility. Specifically, how often do we have those situations where we are aligning same-length sequences (or at least sub-sequences) and what are those situations?

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.

This is an open-ended question that will lead to opinion-based answers, which doesn't fit well with the Stack Exchange format. Could you please change your question to relate to a specific situation or problem? If that can't be done, I recommend instead asking this question on a more discussion-based site like reddit.

Why do you choose the Needleman-Wunsch algorithm?  Are you thinking of using a computer programing/ AI/ML approach to this question and project?  When you mention algorithm, I think more of a Python AI/ML approach.

I see ... the application is if two genes had the same number of indels* at different homologous positions within the gene. This happens but equally one gene might have an indel and another gene might lack it, in which case these genes would be different  lengths.
So just to reiterate your method would find genes with an identical number of indels at different homologous positions in a gene. Its application is protein modelling. This is because indels cause a large change in the entropy of a protein, so even if these two genes had identical sequence (barring the indels obviously) they would be structurally different.
Beyond this however the application is limited because apart from the known changes to protein structure, little is understood about quantifying the function of indels. Certainly in phylogenetics we specifically avoid them, because we struggle to fit them into a known model.

*, indel = insertion deletion event. In a protein they occur in triplets.

My 2 cents: Sequence alignment is a well studied problem in bioinformatics with some great existing methods that do handle uneven length sequences (i.e. indels). With tons of genomic data the trend these days is towards faster methods with more heuristics like pseudo-alignment or k-mer based approaches rather than theoretically optimal but slow and unscalable approaches like Needleman Wunsch. There might be some niche applications on DNA motifs or short peptide analyses but I think these are quite limited.


text:
Fetching uniprot seq using bash command
I am trying to get a protein sequence using this command:
curl https://www.uniprot.org/uniprot/A2Z669.fasta

But I don't get any output:
$ curl https://www.uniprot.org/uniprot/A2Z669.fasta
$

And no file is created.

The URL https://www.uniprot.org/uniprot/A2Z669.fasta redirects to https://rest.uniprot.org/uniprotkb/A2Z669.fasta. You can see this if you open it in your browser: clicking on the first link will take you to the address of the second.
Now, curl doesn't follow redirects by default, you need the -L option for that (see man curl):

-L, --location
(HTTP)  If  the  server  reports  that  the requested page has moved to a different location (indicated with a Location:
header and a 3XX response code), this option makes curl redo the request on the new place. If  used  together  with  -i,
--include or -I, --head, headers from all requested pages are shown.

So you need -L for curl to work here:
$ curl -L https://www.uniprot.org/uniprot/A2Z669.fasta
>sp|A2Z669|CSPLT_ORYSI CASP-like protein 5A2 OS=Oryza sativa subsp. indica OX=39946 GN=OsI_33147 PE=3 SV=1
MRASRPVVHPVEAPPPAALAVAAAAVAVEAGVGAGGGAAAHGGENAQPRGVRMKDPPGAP
GTPGGLGLRLVQAFFAAAALAVMASTDDFPSVSAFCYLVAAAILQCLWSLSLAVVDIYAL
LVKRSLRNPQAVCIFTIGDGITGTLTLGAACASAGITVLIGNDLNICANNHCASFETATA
MAFISWFALAPSCVLNFWSMASR

Or, alternatively, you can use wget which does follow redirects by default:
$ wget "https://www.uniprot.org/uniprot/A2Z669.fasta"
--2023-12-10 16:21:46--  https://www.uniprot.org/uniprot/A2Z669.fasta
Loaded CA certificate '/etc/ssl/certs/ca-certificates.crt'
Resolving www.uniprot.org (www.uniprot.org)... 193.62.193.81
Connecting to www.uniprot.org (www.uniprot.org)|193.62.193.81|:443... connected.
HTTP request sent, awaiting response... 301 Moved Permanently
Location: https://rest.uniprot.org/uniprot/A2Z669.fasta [following]
--2023-12-10 16:21:46--  https://rest.uniprot.org/uniprot/A2Z669.fasta
Resolving rest.uniprot.org (rest.uniprot.org)... 193.62.193.81
Connecting to rest.uniprot.org (rest.uniprot.org)|193.62.193.81|:443... connected.
HTTP request sent, awaiting response... 301 Moved Permanently
Location: https://rest.uniprot.org/uniprotkb/A2Z669.fasta [following]
--2023-12-10 16:21:46--  https://rest.uniprot.org/uniprotkb/A2Z669.fasta
Reusing existing connection to rest.uniprot.org:443.
HTTP request sent, awaiting response... 200 OK
Length: 314 [text/plain]
Saving to: ‘A2Z669.fasta.1’

A2Z669.fasta             100%[=======================================>]     314  --.-KB/s    in 0s

2023-12-10 16:21:47 (8.57 MB/s) - ‘A2Z669.fasta.1’ saved [314/314]

That creates a A2Z669.fasta file in the current directory which will contain the sequence.


text:
Fst calculation from VCF files
I have four vcf files, SNPs_s1.vcf, SNPs_s2.vcf, SNPs_s3.vcf, and SNPs_s4.vcf, which contain information about SNPs. These vcf files were obtained by using the following methods: the initial input files were short-paired reads

I did mapping with minimap2
./minimap2 -ax sr ref.fa read1.fq.gz read2.fq.gz > aln.sam

converted to bam file
samtools view -Sb aln.sam > aln.bam

Sorted the bam file
samtools sort aln.bam > aln_sorted.bam

created an index file
samtools index aln_sorted.bam

created an index file of reference
samtools faidx ref.fa

extracted the raw variants using freebayes
freebayes -f ref.fa aln_sorted.bam > out_var_s1.vcf

extracted SNPs using vcf tools
vcftools --vcf out_var_s1.vcf --remove-indels --recode \
   --recode-INFO-all --out SNPs_s1.vcf

I would like to calculate the fixation index (Fst) between SNPs_s1.vcf and SNPs_s2.vcf, SNPs_s1.vcf and SNPs_s3.vcf etc. I attempted to calculate Fst using bcftools, but I have to provide the id list. I tried the code below for that.
bcftools query -l snps_s1.vcf

The output was "unknown".
How do I calculate Fst using the above vcf files? Your suggestions would be appreciated.

Can you please explain a bit more about Fst, and why you need to calculate it? I'm used to Fst being interpreted as a population statistic, and it doesn't [currently] make sense to me to calculate it for only a few individuals.

@gringer I would like to use Fst to find the genetic diversity between populations. Which method is better for few individuals?

Who is asking you to calculate Fst, and why? If you're doing population genetics, you need populations. A few individuals do not make a population, unless it's kākāpō or kakī, or some other similarly endangered species.

The first issue is easy, the name of your sample in your snps_s1.vcf is set to unknown. You can fix that bit with sed -i '/^#CHROM/s/unknown/sample1/' snps_s1.vcf. But gringer's points are the main issue. What are you really trying to do here with only 4 individuals?


text:
Fastq file with very high per base sequence quality and no box-whisker plot
I have my data from NCBI. I am having difficulty making sense of and correlating my fastq file read 1 and read 2 and the fastQC report. As shown below, I have two paired-end read results each displaying 127 bp in length; however, my results on the fastQC report of per base sequence quality show 91 and it does not show the interquartile range (box-whisker plot). What could be the problem?
Below is an excerpt display of my fastq file and then the fastQC report:
@SRR20708369.1 NS500749:581:HGV7VBGXB:1:11101:17225:1048 length=127
AGCTATCACAGCANGGTACTTCCCAGCGTTACGGGANTGGCCGTACGCCTAACCGCTAACATTACTGCAGGCAACCTACTCATGCACCTAATTGGAAGCGCCACCCTAGCAATATCCACCATTAACC
+
A6AAAE6EA/AAA#EAEAEEEAEEEE/EAEEEEA<E#AAAAEEEE6E<E//EE6EAEEAEAEE/AAA6E//E/EE/A<AEEEE/AEEE//<EEA/EE/E/6<<EE<EEE/<EEEEE//<AE/AAEA/

@SRR20708369.2 NS500749:581:HGV7VBGXB:1:11101:22488:1048 length=127
AGCTATCACGCCANAAGGCTTCCGAGTTAGTATAGGNGCAAAGTGAATACATCTCTGTAACCAGACTCAGATCAAGAAATGATACATTCTGGCCAGGTGTGGTGGCTCATGCCTGTAATCTCAGCAC
+
AAAAAE/EAAAAA#EEEEEEEEEEEEEEEEEEEEE6#AA/AEA/EEEEEEEEAEEEEAEA/EAE/AEAEEA6</6E///EE//6A//AEAEE<E/EE<E<EEEEEEEEEAEAEAAE///EEE/EE6/

:

Cross-posted and answered on biostars: https://www.biostars.org/p/9583849/#9583850

https://trace.ncbi.nlm.nih.gov/Traces/?view=run_browser&acc=SRR20708369&display=metadata is a 10x single-cell sample.
So FastQC plot you are showing us appears to be the R2 read from this sample which is 91 bp (split files, where as the example read you show above seems to be Illumina Index + R1 + R2). This is normal. Illumina index (8 bp), Read 1 is 28 bp (Cellbarcodes + UMI) and R2 is RNA (91 bp) = 127 bp.


text:
What is the best way to acquire protein isoform sequence alignment?

[Update] Thanks @terdon. To clarify my question:
I have a bunch of protein isoforms sequences (produced by the transcripts in the figure) and I want to align protein products (e.g., proteins produced by the first 3 protein-coding transcripts against each other). See there are 3 blocks missing in 5' region in BRCA2-204 compared against BRCA2-206.
This case might be very easy with any MSA softwares such as clustal omega, but there are complicated cases which may trick software into aligning single amino acids instead of being splicing aware. Instead of manually adjusting the alignments, I want to know whether I can directly download, or "translate", this figure into alignment results (e.g., fasta format or clustal format)
This figure should be visible with https://www.ensembl.org/Homo_sapiens/Gene/Summary?db=core;g=ENSG00000139618;r=13:32315086-32400268
in the gene summary tab
[Original] Is there a way to directly download the sequence alignment result from Ensembl? If not, what will be the best way to align exons and introns without parameter-tuning

Please [edit] your question and clarify what you need. You mention protein isoforms in the title, but show transcript sequences in the body. Do you want the alignment of the proteins produced by each transcripts against each other? The alignment of each transcript against a reference genome? The alignment of the transcripts against each other? Ideally, give us an example we can look up, a link that we can click on to see what you are using (is that human BRCA2 or something else?) and then explain what you would like to get from that.

@terdon Thanks!

Interesting question (upvoted). This might be naive, but if these are transcripts surely they have already been spliced? That aside from eukaryotic genomics its easy, but must have been solved years ago, I suspect a de novo pipeline is simply "reinventing the wheel".

@M__♦ I am not sure I followed. You mean aligning these spliced exons is easy?

I've undeleted the question ... it's a perfectly good question - it has an answer by an expert which is upvoted. My justification is that I had trouble logging in over the last 24 hours - I had to escalate to SE login HQ to get in, hence I couldn't take part in the discussion before the question was deleted. @terdon is an expert on human genetics and I'll defer to their answer.

Please don't deface your question. Once a question has been posted and answered, you have given SE the right to use the content under CC-BY-SA. Even if your question is no longer useful to you, or if you have your answer, that doesn't mean that other users won't also find it useful so deleting useful content from the site isn't allowed.

Here is a manual (or semi-manual, at best) approach. If you need to do this for very many sequences it won't be practical, but if dealing with one or two genes, it might be enough. The idea is to:

Get the protein coding transcripts
Get the proteins they code for
Align

In the page you gave, there is a "Show transcript table" button near the top:

If you click on that, you get the list of transcripts and their corresponding proteins:

You can then click on the "Show/hide columns" and select "Translation ID" so you can get the Ensembl protein corresponding to each transcript:

If you now deselct everything and only leave "Translation ID" selected:

You are left with a nice list of Ensembl proteins. Copy/paste those into a text file:
ENSP00000369497.3
ENSP00000499438.2
ENSP00000505508.1
ENSP00000439902.1
ENSP00000514856.2
ENSP00000434898.2
ENSP00000506251.1
ENSP00000499570.2
ENSP00000433168.2
ENSP00000499256.2
ENSP00000514855.1

Now, go to BioMart and retrieve the sequences. Select "Ensembl Genes" as the database and then "Human Genes":

Then, click on "Filters" (1) and then "Gene" (2), select "Input external references ID list [Max 500 advised]" (3), set the ID to "Protein stable ID(s) with version" (4) and paste the protein list (5):

Next, go to "Attributes" (1) and select "Sequences" (2):

Now, finally, click on "Results" (1) and you will get the sequences in fasta format (2):

I'm afraid they are too long to fit in this answer, but I posted them here: https://pastebin.com/NZycAmgD.
You can now paste these sequences into your favorite aligner and get the alignment you were after. Here, I used the online clustalΩ (the image only shows the first few lines):

You can find the whole alignment in clustal format (text) here.

Thanks @terdon. This is exactly my question. In your clustal omega result screenshot (last figure of your answer) there are a bunch of misalignments for ENST00000700201 compared with others (starting from loc ~105 to the end). I assume the software is forcing different exons aligned together to reduce number of gaps. I mean those gaps should be inserted before loc ~105 on ENST00000700201, not after. But somehow the software prefer before (due to parameters, misalignments outperformed gaps). I can manually adjust it. Easy for one gene, but not for hundreds of genes. That's why I want to downloa

although ENST00000700201 is NMD, not a perfect example. But the question stays for other protein coding transcripts.

@JenMarylin I don't quite understand what you mean. What exons? We're aligning simple protein sequences here, the aligner has no knowledge of the exonic structure of the transcripts that were translated to make these proteins. And this is a multiple sequence alignment, it is positioning the gaps so as to maximize the score of the entire alignment. What is it you want to download? Are you looking for an alignment of the translated exons instead of the actual protein product?

This is a paper that tries to address this problem
"Splice-Aware Multiple Sequence Alignment of Protein Isoforms"
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6508070/
I have not evaluated it, but it could be a point of reference
note that there is also a sequel to this paper from 2023 https://www.ncbi.nlm.nih.gov/pmc/articles/PMC10166558/