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

    }

}


text:
I am trying to find correlation between two ChIP-seq datasets but the Pearson and Spearman test results do not correspond with the peak overlap
I am trying to compare two ChIP-seq datasets for histone modifications, however, the Spearman and Pearson correlations between two particular datasets does not correspond to their peak overlap (peak overlap is very low, while both pearson and spearman tests are showing strong correlation). I am performing this analysis on Usegalaxy, first using multibamSummary on both bam files followed by plotCorrelation. What could be the error here?

Can you please show any plots or data that have helped you come to the conclusion that there is a problem?

I can readily imagine a situation where the bulk data agree very well in terms of the overall genomic profile, but the statistical process by which you pick peaks yields different results. Knowing more about the samples would be handy- are they technical replicates, biological replicates, experimental arms?

I agree with this, the correlation could be supported by the breadth of values across the data set. The actual r2-value is key: if it's 0.5 to 0.7 I wouldn't consider there to be a discrepancy. If it's >0.8 thats puzzling. At >0.9  there's definite selective sampling occurring because thats weird. The r alone (Spearman's) will be higher in most cases so I would use r2 as the guide. Summary a correlation can be stabilised by agreement between low and high values, so the exact correlation score is needed.

Correlations like this between ChIP-seq datasets is common.
If you're calculating the correlation between BAM files across the entire genome, those are largely going to be the same (hence the need for input controls in ChIP-seq for performing the peak calling in the first place).
Similarly, many ChIP-seq experiments only have 30 - 60% of reads within their peaks (FRiP), which means that 40 - 70% of reads don't fall within peak regions.
This is going to heavily skew any correlation calculations that don't exclusively focus on the peaks themselves.
To more clearly compare ChIP-seq data against each other, global metrics like FRiP can help you diagnose measurement/sequencing issues.
Looking at the correlation across ChIP-seq signal within peaks may be a more informative calculation than looking at all reads across the genome.
Similarly, principal component analysis on signal within peaks can be helpful here.
DiffBind has some helper functions to perform these types of calculations that you can look at.

source: Stack Exchange

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

sort: votes


text:
variant visualisation of VCF files for 30+ samples
I am relatively new to genome analysis and would like to compare 30+ vcf files, 1 sample per file. After filtering there are about ~20 variants per vcf file.
I would like to evaluate to mutational landscape between samples and for each sample as well.
One option that I found is a mutation heatmap where each row is a sample and each column a position in the genome. Is this the preferred way of displaying this type of information and how would I include gene information?
At this point I am interested in getting an overview first of the affected positions and then dig deeper using e.g. zygosity.
Thank you very much for helping out!
2nd edit
I would also like to evaluate if there is a bias towards certain genes. For example are certain oncogenes mutated more frequently in this cohort or is the accumulation random and unique within each sample?
Best,
M

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.

What exactly do you want to compare? Simply whether a variant was called in the different samples? Would you want to take zygosity into account (is it a match if the same variant is found heterozygous in one sample and homozygoyus in another)? Do you want to compare the genes which had variants found in them and how they segregate in your samples? What exactly does "visualize a variant" mean?

Thank you very much for the comments and I updated the question. Currently I am searching for a good way to display the position of a variant in the genome for each sample and include e.g. gene annotation. @terdon Then I would like to compare the affected genes, for example are the mutations biased towards oncogenes.

Thanks for the extra context. I'll reiterate that you cannot assume that the lack of a variant in a PVCF file implies that the sequence is homozygous reference at that locus. That locus could be poorly sequenced and therefore be a "no call". You need to "joint call" the original GVCF files to be confident.
OK, but your question: I don't work in cancer, so I'm unfamiliar with what you're trying to do. It sounds like you have ~60 genotypes. Both R and Python have libraries for reading VCFs and visualization. The R library ggplot has gemo_tile which can create heat maps.

It's a bit hard for me to advise about your second edit because it involves questions of science (I'm a software engineer). Do you have a particular statistical distribution of variation within genes that you think is normal? Again, these questions might be obvious to someone working in cancer, but I'm unfamiliar with standard practice there. Perhaps a two-sample Student's t-test is appropriate? With just 60 genotypes, I'd just manually copy each one from the VCF and into an R script and use the R Student's t-Test function: https://stat.ethz.ch/R-manual/R-devel/library/stats/html/t.test.html .

Single sample PVCF files are generally not combinable in the way you describe. PVCF stands for "project" VCF. As far as I know, PVCF isn't defined anywhere but, it refers to a multi-sample VCF that does not contain reference blocks. It only contains lines for which at least one sample has at least one alternate allele. We distinguish this from a GVCF which is a single-sample VCF that does contain reference blocks. We can't combine two PVCFs in a principled way because some variant might only appear in one PVCF. We lack the reference block information necessary to decide if a sample in the other VCF is 0/0 or "no call".
Do all the VCF files have the same set of variants? If all the VCF files have the same set of variants, you might try using Hail (disclaimer: I'm a Hail maintainer):
import hail as hl

vcf_files = ['sample1.vcf.bgz', ...]
mts = [hl.import_vcf(f) for f in vcf_files]
mt = hl.MatrixTable.union_cols(mts)  # assumes each file has the same set of variants
mt.write('combined.mt')  # save to Hail's native format so the following code is fast

mt = hl.read_matrix_table('combined.mt')
mt.show(n_rows=5, n_cols=10)  # see a table of the first 5 rows and 10 cols

mt = hl.variant_qc(mt)  # adds a bunch of variant QC annotations
mt.describe()  # inspect available annotations

You mentioned "genomic distribution". Are you referring to the allele frequency in this sample? You can see a genome-wide scatter plot of reference allele frequency like this:
from hail.ggplot import *
fig = ggplot(mt) + geom_point(aes(x=mt.locus, y=mt.variant_qc.AF[0]))
fig.show()

Thank you for your answer. Each vcf file is derived from a different sample and I expect minimal overlap between the called mutations. I updated my post, but in brief, I would like to summarise the mutational burden between samples and evaluate if there is a bias towards certain genes. for example are certain oncogenes mutated more frequently in this cohort or is the accumulation more random?

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


text:
calculating mutation frequencies for every gene
I have a dataset for mutation data and I want to calculate mutation frequencies across all genes
df (This is only the small subset of data)
Gene name   Sample id   MUTATION_ID Mutation Description
ARID1B  2719660 171258500   Substitution - Missense
ARID1B  2719659
ARID1B  2719661 171258501   Substitution - Missense
ARID1B  2719662
ARID1B  2719663 171258501   Substitution - Nonsense
CD58    2878555 110346783   Substitution - Nonsense
CD58    2877956
CD58    2878557
CD58    2877958 110346784   Substitution - Nonsense
CD58    2878559 110346785   Substitution - Nonsense
CD58    2877960
MRE11   2861617 123320443   Substitution - coding silent
MRE12   2861617 123320444   Substitution - coding
MRE13   2861617
MRE14   2861617 123320445   Substitution - coding silent
MRE15   2861617
MRE16   2861617 123320446   Substitution - coding

The formula for calculating the mutation is
Positives ÷ (Positives + Negatives) x 100
where,
Positives  = No of samples where MUTATION_ID is present
Negative = No of samples MUTATION_ID for the sample
I want to calculate mutation frequency for every gene in the column_1:Gene name with python script
I tried the following code
df = df.groupby("Gene name").count()
Positives = df["MUTATION_ID"]
Negatives = df["Sample id"] - df["MUTATION_ID"]
df['Mutation_Frequency'] = Positives / (Positives + Negatives) * 100

Are your fields tab separated? They probably have to be since you have spaces in the description field. And what are those extra -? Is that a separate field with no header?

The code is fine. The only issue is you'd put \s characters in-between the column header names. I've place an additional "mutation_frequency" column because thats how I'd do it.
df = pd.read_csv('/pathtodir/test/bioinfo.csv',sep="\s",engine='python')
df = df.groupby("Genename").count()
df['Mutation_Frequency1'] = ((df["MUTATION_ID"] / (df["Sampleid"] + (df["Sampleid"] - df["MUTATION_ID"]))) * 100).round(2)
df['Mutation_Frequency2'] = (df["MUTATION_ID"] / df["Sampleid"])*100
df.drop(['Description', 'Mutation'], axis=1, inplace=True)
print(df)

-
Sampleid
MUTATION_ID
Mutation_Frequency1
Mutation_Frequency2

ARID1B
5
3
42.86
60.0

CD58
6
3
33.33
50.0

MRE11
1
1
100.00
100.0

MRE12
1
1
100.00
100.0

MRE13
1
0
0.00
0.0

MRE14
1
1
100.00
100.0

MRE15
1
0
0.00
0.0

MRE16
1
1
100.00
100.0


text:
Inconsistent replicate numbers in RNA-seq
I have 3 groups for the RNA-seq analysis (Control, treatment A and treatment B). There are 2 replicates for control and treatment A and 3 replicates for treatment B (lost 2 replicates due to a mistake in sample handling). May I ask how will the inconsistent replicate numbers impact the results (e.g. DEGs)?

It will not do any harm other than a reduction in statistical power compared to more samples per group. Common analysis tools such as DESeq2, edgeR and limma will work on unequal sample size just fine.


text:
How to align multiple contigs? To generate full length genome
I have MAG geneome sequence released in NCBI. In the released sequence there are 4915 contigs present. I want to get the full length genome sequence as FASTA file to proceed with further works like circular genome viewing and other molecular typing. Kindly let me know how to proceed in this regard. I am not much aware of bioinformatics so need even the basic information in this regard.
Thanks

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

Please explain what genome you are working with in particular and how you have assembled this to date. Circular genome viewing ... bacteria?

Just to note @Laura is an expert on this subject - if its bacteria. This post might help https://bioinformatics.stackexchange.com/questions/20336/improving-prokaryotic-assembly-with-other-contig-scaffold-level-data

Simple answer is, it might not be possible to do it.
The large number of contigs suggests to me that a short-read sequencing method (such as Illumina sequencing) was used. Alone, short-read sequencing technologies cannot resolve repetitive sequences that are longer than the length of the reads. This is why you have a large number of contigs.*
*You might already have done this, but if you did not, it might be worth checking that your MAG is not actually a mixture of different genomes. This can also somewhat increase the contig number and, more importantly, can skew your interpretation if you assume it's only one organism. There are various ways to test this, but the simplest I know is the program CheckM. To run it you will need a way to use Unix/Linux (could be on a virtual machine). If you are not familiar with Unix, bash and conda, I would heartily suggest software carpentry lessons such as The Unix Shell (https://swcarpentry.github.io/shell-novice/) and possibly this conda lesson although it is in beta and I haven't tested it https://carpentries-incubator.github.io/introduction-to-conda-for-data-scientists/
I'm not familiar with the term 'molecular typing', but it looks like you're interested in what organisms / genes there are in your sample? You probably might not need a full genome for that.
Small probability of helping - but it might also be good to have a look at the sizes of your contigs. You might luck out and have a few large contigs (hundred thousand - million basepairs) which could be compared to existing database genomes.


text:
How to extract a list of hydrogen bonds from a PDB?
I am analyzing the 7VF2 structure, from the RCSB-database. I noticed that some residues belonging to one of its chains interact with the other chains that are present within the assembly. My question is if there is a method to extract a list of interacting residues that are present in a PDB.
EDIT:

Just as an example, we have Lys 1616 belonging to ZC3H13 chain that interacts with Glu 621 of the protein virilizer homolog (the first chain)
Mol* 3D visualizer can be used for major information in the 3D View in the link page I provided

do you need alpha chain hydrogen interactions or side chains hydrogen interactions ? I believe most of the alpha helics and beta sheets emerge from apha carbon hydrogen bonds so thet'll be a huge number of interaction to be shown, but could be wrong

https://files.rcsb.org/pub/pdb/validation_reports/vf/7vf2/7vf2_full_validation.pdf : Atom-1 Atom-2 Interatomic distance (Å) Clash overlap (Å) : 2:B:1616:LYS:HA 2:B:1616:LYS:HD2 1.71 0.42 there is something wrong with k1616 as per validation report states ... see atoms nomenclature in Lys http://ligand-expo.rcsb.org/reports/L/LYS/LYS_D3L3.gif

DSSP gives you hydrogen bonds, not sure how it works for complexes though.

You can do that with GROMACS. You probably need to check the manual for this, but roughly these are the steps:
gmx pdb2gmx -f input.pdb -o output.gro -water spc
gmx solvate -cp output.gro -cs spc216.gro -o output_solv.gro -p output.top
gmx pdb2gmx -f output_solv.gro -o output_processed.gro -p output.top -ignh
gmx make_ndx -f output_processed.gro -o index.ndx
gmx hbond -f output_processed.gro -n index.ndx -s output.tpr -num hbond.xvg -hbn hbond.ndx

First, you convert the PDB file into a GROMACS topology file using the pdb2gmx command. Then you can solvate the system with water molecules and add hydrogen atoms. Then you cCreate an index file containing the groups of atoms you want to analyze. Then you can use the hbond command to calculate the hydrogen bonds.
I have not used GROMACS in a while. Hope that provides a starting point. With GROMACS you can even run a simulation and monitor the h-bonds over time.
You actually might need to run an energy minimziation first, as the PDB might not in an optimal state for the h-bonds detection.

Great comment! I will keep this question open so that more possible answers can be added. Meanwhile, thank you for the suggestion! I will try and let you know!


text:
Which basecaller for nanopore is the best to produce event tables with information about the block size/move table?
I previously used albacore version 2.3.1 to make initial move tables, but then I re-squiggle using Tombo version 1.5.1 to fix the errors.
Example of move table produced by albacore, in which all the moves have the same size:

When it is re-squiggled by Tombo (which produces variable length moves/blocks):

I am starting a new project, and I can repeat the same pipeline, however, my concern is that albacore version 2.3.1 is too old now and I was wondering to ask if there are better pipelines/software for base-calling that produces initial move tables so that I can build on top of.

ONT's newest basecaller, Dorado, can do this for SAM output with the --emit-moves argument. It was added in this commit.
Dorado can be downloaded from here:
https://github.com/nanoporetech/dorado
Here's an example command line:
$ ~/install/dorado/dorado-0.2.1-linux-x64/bin/dorado basecaller --emit-moves \
  ~/install/dorado/dorado-0.2.1-linux-x64/dna_r10.4.1_e8.2_260bps_sup@v4.1.0 DAE_pod5 | \
  samtools view -b - > called_DAE_gDNA_reads_all.bam

And a few lines of example output:
@HD     VN:1.6  SO:unknown
@PG     ID:basecaller   PN:dorado       VN:0.2.1+c70423e        CL:dorado basecaller --emit-moves /home/minion/install/dorado/dorado-0.2.1-linux-x64/dna_r10.4.1_e8.2_260bps_sup@v4.1.0 DAE_pod5
@PG     ID:samtools     PN:samtools     PP:basecaller   VN:1.10 CL:samtools view -b -
@PG     ID:samtools.1   PN:samtools     PP:samtools     VN:1.10 CL:samtools view -h called_DAE_gDNA_reads_all.bam
@RG     ID:460fc479-3b30-4ed5-801f-a82bbe556101_dna_r10.4.1_e8.2_260bps_sup@v4.1.0      PU:PAK83314     PM:     DT:2023-02-08T10:22:53.256+00:00        PL:ONT  DS:basecall_model=dna_r10.4.1_e8.2_260bps_sup@v4.1.0 runid=460fc479-3b30-4ed5-801f-a82bbe556101 LB:DAE_gDNA     SM:DAE_gDNA
@RG     ID:227379b1-c5fc-465d-9e24-eb771133f11e_dna_r10.4.1_e8.2_260bps_sup@v4.1.0      PU:PAK83314     PM:     DT:2023-02-08T20:28:34.34+00:00 PL:ONT  DS:basecall_model=dna_r10.4.1_e8.2_260bps_sup@v4.1.0 runid=227379b1-c5fc-465d-9e24-eb771133f11e LB:no_sample    SM:no_sample
003d84eb-ce0f-47ea-9c7c-fe939fe8702e    4       *       0       0       *       *       0       110     GTGCCTCCTTGCTTTGGTTCAGTTCTATCCAGCCTTTAGATATCCCAGCCCTCACAGGGTTTGTGTGACCTTGTTGCATTATCATGCGAAACACTTTCGCGTTTTTAGTG  %%%%%(&%%&%%'((''))*.-,-''&&&'&&&(%$$$%&&&'%%%%&(*))'&&&&)-05>6;>>777887++'&&&&&&%%%&*.....0820001:;;;<<87''&&  qs:i:7  du:f:0.526      ns:i:2104       ts:i:10 mx:i:4  ch:i:1037       st:Z:2023-02-09T07:07:00.271+00:00      rn:i:96929      f5:Z:PAK83314_fail_unclassified_227379b1_ea81734c_853.pod5      sm:f:98.674     sd:f:22.395     sv:Z:quantile   RG:Z:227379b1-c5fc-465d-9e24-eb771133f11e_dna_r10.4.1_e8.2_260bps_sup@v4.1.0    mv:B:c,5,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,1,0,0,1,1,0,0,1,0,0,1,0,0,0,0,1,0,1,1,0,1,0,1,0,1,0,1,0,0,1,0,1,1,0,0,1,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,1,0,1,0,1,0,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,1,0,0,1,0,0,0,0,0,1,0,1,0,1,1,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,0,1,0,0,0,1,0,0,1,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,1,0,0,1,0,1,0,1,0,0,1,0,1,0,0,0,1,0,0,0,1,0,1,0,1,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,1,0,1,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,1,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,1,0,0,1,0,1,1,0,0,1,0,1,0,0,1,0,1,0,1,0,0,0,0,1,0,0,0,1,0,0,1,0,0,1,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,1,0,1,0,0,1,0,1,0,1,0,1,0,0,1,0,0,1,0,1,0,0,0,1,0,1,0,1
002e48ca-84bf-4419-8370-943d1ca3ac64    4       *       0       0       *       *       0       107     TATGTCTAACTTACTTCGTTCAGTTGTGAAAACTGGGTGTTTAACCTTTGATTCTATCGTGTTTCCCTAGTTTTCGCATTTGACGTGAAACGCTTTCACATTTTGGT     &&'''&%%%%((((*+)***())((')&%%%%$%%%'--,,,,-6996/,+-5ABBBCCFBCE>88:1764424/--,,--&%$$%*08??;;;;<B<:((('()//     qs:i:9  du:f:0.92975    ns:i:3719       ts:i:10 mx:i:4  ch:i:1058       st:Z:2023-02-09T07:07:08.905+00:00      rn:i:62167      f5:Z:PAK83314_fail_unclassified_227379b1_ea81734c_853.pod5      sm:f:85.094     sd:f:23.944     sv:Z:quantile   RG:Z:227379b1-c5fc-465d-9e24-eb771133f11e_dna_r10.4.1_e8.2_260bps_sup@v4.1.0    mv:B:c,5,1,0,0,0,0,0,0,1,0,0,0,0,0,1,1,0,0,0,0,0,1,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,1,1,1,0,1,0,1,0,0,1,0,0,1,0,1,1,0,0,1,0,1,0,0,0,0,1,0,0,1,0,1,0,0,1,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,1,0,1,0,0,1,0,1,0,1,1,0,1,0,0,0,0,0,1,0,0,1,0,0,1,0,1,0,1,0,0,1,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,1,0,0,0,0,0,0,1,0,1,0,1,0,0,1,0,1,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,1,0,1,0,0,1,0,1,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,1,0,1,0,1,0,1,0,1,0,0,1,0,1,0,0,0,1,0,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,0,1,0,0,1,0,0,1,0,0,0,0,1,0,1,0,0,1,0,0,0,1,0,0
0138d83a-da9f-4468-bafc-b974c4e9eece    4       *       0       0       *       *       0       234     ATGTTTATGTGTTCTTCGTTAGTAAGCCACCGATATTGACTCCAGGGTTTGTGTAACCTTGTTTTCGCATTTATCGTGAAACGCTTTCGCGTTTTTCGTCGCTTCAAATAGTTTTAATTCAGTTGTCTAAAAATGTTAACATATGGGTATTATTTCTCTAGAACATTTTAGTTGTGTAATTTCTACCATTTAAATGGGGTGAATAAATAACCAAATCATTGAACAACTGGAATG      &&&')'&%$$##$&&''+,-)''&'')**'&%&''))&&&)'((*+8<AEEFGHJIJIJIICA@;843344:254,,,,CEE>====7<<<>>@=?+***'&&''(''&%%%%$%%&)1//..-+,+*'&'(+2+++45<ACCC66677@75556+((((55>@AEFHFHGFEHHJHHFHIKJGFFFKLHZEGEGGIDFDDCDDERKIEGDG>=>5332/'&&&&)**,2242.      qs:i:10 du:f:0.69975    ns:i:2799       ts:i:10 mx:i:4  ch:i:1933       st:Z:2023-02-09T07:07:01.924+00:00      rn:i:77788      f5:Z:PAK83314_fail_unclassified_227379b1_ea81734c_853.pod5      sm:f:97.397     sd:f:27.276     sv:Z:quantile   RG:Z:227379b1c5fc-465d-9e24-eb771133f11e_dna_r10.4.1_e8.2_260bps_sup@v4.1.0    mv:B:c,5,1,0,0,0,1,0,1,0,1,0,1,0,0,1,1,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,1,1,0,0,1,0,0,0,1,0,0,1,0,1,0,0,1,0,1,0,0,1,0,0,0,0,0,0,1,0,1,0,1,0,1,0,1,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0,0,1,0,1,1,0,0,1,1,0,0,0,1,0,0,0,1,0,0,1,0,0,1,0,1,0,0,0,0,1,1,0,0,1,0,1,0,0,1,0,1,0,0,0,1,0,1,0,1,0,0,1,0,1,0,0,1,0,0,1,0,1,0,0,1,0,0,1,0,0,1,0,0,0,1,0,1,0,1,0,1,1,0,0,1,0,1,0,1,0,1,0,0,1,0,1,0,0,1,0,1,0,1,0,0,0,0,0,0,0,1,1,0,0,0,1,0,1,0,1,0,1,0,0,0,1,0,1,0,1,0,1,0,1,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,1,1,0,1,1,0,1,0,1,1,0,1,0,1,0,0,0,0,0,0,1,0,1,0,1,0,1,0,1,1,1,0,1,0,1,0,1,1,1,1,0,1,0,1,0,1,0,1,0,1,0,1,0,0,1,1,0,1,1,0,1,0,1,1,1,1,0,1,0,1,1,1,0,0,1,1,0,0,0,0,1,0,1,0,1,0,0,1,0,0,0,0,0,1,0,1,0,1,0,1,0,1,0,1,1,0,1,0,1,1,1,1,0,0,1,0,1,1,0,1,0,1,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,0,0,0,1,0,1,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,0,1,0,1,0,0,1,0,1,0,0,1,0,0,1,0,1,0,1,0,1,0,1,0,1,0,1,1,0,1,0,1,0,1,0,1,0,0,1,0,0,1,0,0,1,0,1,1,0,1,0,1,0,1,1,0,1,1,0,1,1,0,1,0,1,0,0,1,0,0,1,1,0,0,1,0,0,1,0,0,1,1,1,1,0,0,1,0,1,0,1,0,1,0,1,0,1,1,0,1

Based on the commit, the move table is attached to the mv tag for each read, indicating an array of 8-bit signed integers (represented as comma-separated values in the SAM output). I'd guess (based on how the calling has worked in the past) that these represent strides comprising a fixed number of consecutive samples, so (as before) it'll need correction to convert the stride width to variable width events.
For the model I'm using here, my guess from the model toml file is that the stride width is 5:
$ grep -B 8 stride ~/install/dorado/dorado-0.2.1-linux-x64/dna_r10.4.1_e8.2_260bps_sup@v4.1.0/config.toml
[encoder]
type = "serial"
[[encoder.sublayers]]
type = "convolution"
insize = 1
size = 16
bias = true
winlen = 5
stride = 1
--
max = 3.5

[[encoder.sublayers]]
type = "convolution"
insize = 16
size = 16
bias = true
winlen = 5
stride = 1
--
max = 3.5

[[encoder.sublayers]]
type = "convolution"
insize = 16
size = 1024
bias = true
winlen = 19
stride = 5
```


text:
Can you recommend me some programs or R packages for species distribution modelling?
I am working with a whole genus of land snails, trying to predict its distribution in the coming decades and how it might get affected by climate change. I have already worked with MaxENT and the R package Biomod2 (I have tried out all the important algorithms within it as well). I was wondering if any of you can recommend another option thats reliable.

Anything with the negative binomial distribution will always work.

It might be helpful to know what ecological data you are working from.

Not sure if that's what you're asking but I'm working with presence only occurrence data, climate data from worldclim and I also possess phylogenetic data but I have yet to use them in such work. Also @M__ id like it if you were a little more specific and maybe give some reasoning behind your response.

Spatial species distribution always follow a negative binomial distribution at any point of most other factors e.g. climate. When a species enters a human monoculture type environment - then deviation from the negative binomial can occur within that monoculture. The clumping will be approximated via trigonometry: 'cause the clumping can be approximated via a circle, very clumped - small circle, more dispersed bigger circle.
The parameters of the negative binomial are likely to change in relationship to climate and thats what I'd model if I was doing this work.  The 'clumping' (which is what the negative binomial models) is likely to tightened as global warming increases. In practical terms the snails start drying out so they have to clump in the remaining pockets of moisture (this might be more true of slugs than snails ... at a guess). Alot of invertebrates that cannot or don't really fly - are prone to desiccation.
Personally, I go to empirical stats rather than use a package. However, the packages you are using will follow these two principles. So you can do it from first principles. Its fairly simple ecology its not like phylogenetics that have been through decades of development.
If you've DNA from the snails then the phylogenetic Bayesian spatial models around Beast provide a powerful basis for the investigation. That would be very strongly recommended if you've DNA sequence from a good proportion of the snails.


text:
Sanger Sequencing Knitting Error
I am doing a project where I am reproducing the analysis from the article "sangeranalyseR: Simple and Interactive Processing of Sanger Sequencing Data in R". Below is the example chunk for Sanger Alignment. Below are the package and datasets.
library(sangeranalyseR)
data("sangerAlignmentData")
data("sangerContigData")
data("sangerReadFData")

Below is the code for Sanger Alignment, similar to the code for Sanger Contig and Sanger Reading.
rawDataDir <- system.file("extdata", package = "sangeranalyseR")
fastaFN <- file.path(rawDataDir, "fasta", "SangerAlignment", "Sanger_all_reads.fa")
my_sangerAlignmentFa <- new("SangerAlignment",
                            inputSource          = "FASTA",
                            processMethod        = "REGEX",
                            FASTA_File           = fastaFN,
                            REGEX_SuffixForward  = "_[0-9]*_F$",
                        REGEX_SuffixReverse  = "_[0-9]*_R$",
                            refAminoAcidSeq      = "SRQWLFSTNHKDIGTLYFIFGAWAGMVGTSLSILIRAELGHPGALIGDDQIYNVIVTAHAFIMIFFMVMPIMIGGFGNWLVPLMLGAPDMAFPRMNNMSFWLLPPALSLLLVSSMVENGAGTGWTVYPPLSAGIAHGGASVDLAIFSLHLAGISSILGAVNFITTVINMRSTGISLDRMPLFVWSVVITALLLLLSLPVLAGAITMLLTDRNLNTSFFDPAGGGDPILYQHLFWFFGHPEVYILILPGFGMISHIISQESGKKETFGSLGMIYAMLAIGLLGFIVWAHHMFTVGMDVDTRAYFTSATMIIAVPTGIKIFSWLATLHGTQLSYSPAILWALGFVFLFTVGGLTGVVLANSSVDIILHDTYYVVAHFHYVLSMGAVFAIMAGFIHWYPLFTGLTLNNKWLKSHFIIMFIGVNLTFFPQHFLGLAGMPRRYSDYPDAYTTWNIVSTIGSTISLLGILFFFFIIWESLVSQRQVIYPIQLNSSIEWYQNTPPAEHSYSELPLLTN",
                            minReadsNum           = 2,
                            minReadLength         = 20,
                            minFractionCall       = 0.5,
                            maxFractionLost       = 0.5,
                            geneticCode           = GENETIC_CODE,
                            acceptStopCodons      = TRUE,
                            readingFrame          = 1,
                            processorsNum         = 1)

launchApp(my_sangerAlignmentFa)

writeFasta(my_sangerAlignmentFa,
           outputDir         = tempdir(),
           compress          = FALSE,
           compression_level = NA,
           selection         = "all")

generateReport(my_sangerAlignmentFa,
               outputDir           = tempdir(),
               includeSangerRead   = TRUE,
               includeSangerContig = TRUE)

There are two problems that occurred for me. One is that the launchapp() function created an error for me stating
INFO [2023-03-01 14:25:52] Your input is 'SangerAlignment' S4 instance
INFO [2023-03-01 14:25:52] SangerAlignment with 'FASTA' inputSource
cannot run Shiny app (You don't need to do trimming or base calling)
NULL

How can I fix this error?
Another error that occurred within knitting. I ran the same code without the launchapp() function and an error that I never seen before said this
Quitting from lines NA-219 (SangerRead_Report_fasta.Rmd)  Error in parse_block(g[-1], g[1], params.src, markdown_mode) :
Duplicate chunk label 'unnamed-chunk-1', which has been used for the chunk:
library(sangeranalyseR) library(rentrez) Calls: <Anonymous> ...
process_file -> split_file -> lapply -> FUN -> parse_block

Execution halted

Is it my RStudio or is the code itself. Is there anything I am missing?
sessionInfo(sangeranalysesR)
R version 4.2.2 (2022-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 22621)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.utf8
[2] LC_CTYPE=English_United States.utf8
[3] LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.utf8

attached base packages:
character(0)

other attached packages:
[1] sangeranalyseR_1.8.0

loaded via a namespace (and not attached):
  [1] utils_4.2.2            nlme_3.1-161
  [3] fs_1.6.1               bitops_1.0-7
  [5] bit64_4.0.5            httr_1.4.4
  [7] GenomeInfoDb_1.34.6    tools_4.2.2
  [9] bslib_0.4.2            utf8_1.2.2
 [11] R6_2.5.1               DT_0.27
 [13] DBI_1.1.3              lazyeval_0.2.2
 [15] BiocGenerics_0.44.0    colorspace_2.1-0
 [17] rhdf5filters_1.10.0    ade4_1.7-20
 [19] withr_2.5.0            tidyselect_1.2.0
 [21] gridExtra_2.3          phangorn_2.11.1
 [23] bit_4.0.5              compiler_4.2.2
 [25] cli_3.5.0              datasets_4.2.2
 [27] shinyjs_2.1.0          plotly_4.10.1
 [29] ggdendro_0.1.23        base_4.2.2
 [31] sass_0.4.5             scales_1.2.1
 [33] quadprog_1.5-8         stringr_1.5.0
 [35] digest_0.6.31          rmarkdown_2.20
 [37] sangerseqR_1.34.0      XVector_0.38.0
 [39] pkgconfig_2.0.3        htmltools_0.5.4
 [41] sessioninfo_1.2.2      fastmap_1.1.0
 [43] excelR_0.4.0           grDevices_4.2.2
 [45] htmlwidgets_1.6.1      rlang_1.0.6
 [47] rstudioapi_0.14        RSQLite_2.2.20
 [49] shiny_1.7.4            jquerylib_0.1.4
 [51] generics_0.1.3         jsonlite_1.8.4
 [53] BiocParallel_1.32.5    zip_2.2.2
 [55] dplyr_1.1.0            RCurl_1.98-1.10
 [57] magrittr_2.0.3         GenomeInfoDbData_1.2.9
 [59] biomformat_1.26.0      Matrix_1.5-3
 [61] Rhdf5lib_1.20.0        Rcpp_1.0.10
 [63] munsell_0.5.0          S4Vectors_0.36.1
 [65] fansi_1.0.4            logger_0.2.2
 [67] DECIPHER_2.26.0        ape_5.6-2
 [69] shinycssloaders_1.0.0  lifecycle_1.0.3
 [71] yaml_2.3.7             stringi_1.7.12
 [73] MASS_7.3-58.2          zlibbioc_1.44.0
 [75] rhdf5_2.42.0           plyr_1.8.8
 [77] grid_4.2.2             blob_1.2.3
 [79] parallel_4.2.2         promises_1.2.0.1
 [81] shinydashboard_0.7.2   crayon_1.5.2
 [83] methods_4.2.2          lattice_0.20-45
 [85] Biostrings_2.66.0      zeallot_0.1.0
 [87] knitr_1.42             pillar_1.8.1
 [89] igraph_1.3.5           seqinr_4.2-23
 [91] reshape2_1.4.4         codetools_0.2-19
 [93] stats4_4.2.2           fastmatch_1.1-3
 [95] reprex_2.0.2           BiocVersion_3.16.0
 [97] glue_1.6.2             evaluate_0.20
 [99] BiocManager_1.30.20    data.table_1.14.6
[101] vctrs_0.5.2            httpuv_1.6.8
[103] graphics_4.2.2         gtable_0.3.1
[105] purrr_1.0.1            tidyr_1.3.0
[107] cachem_1.0.6           ggplot2_3.4.0
[109] openxlsx_4.2.5.1       xfun_0.36
[111] mime_0.12              xtable_1.8-4
[113] later_1.3.0            viridisLite_0.4.1
[115] tibble_3.1.8           stats_4.2.2
[117] memoise_2.0.1          IRanges_2.32.0
[119] shinyWidgets_0.7.6     ellipsis_0.3.2
[121] BiocStyle_2.26.0
```

Your code runs without error on my laptop (macOS 13.2.1; R v4.2.2; Rstudio 2022.07.2), so I have some further questions that might help with troubleshooting (too many to include in a comment):

Does the example from the docs run without errors/warnings on your system?
Did you install the sangeranalyseR package via Bioconductor (per the docs)?
Do you have the most up-to-date Bioconducter installed?
Have you tried running the example code in a 'fresh' session?
Are you able to please run the command sessionInfo() and edit your question to include the output?
If you run your example via the reprex package, are there any differences between your output and my output?

The code from your question run on my laptop:
# BiocManager::install("sangeranalyseR")
library(sangeranalyseR)
#> Loading required package: stringr
#> Loading required package: ape
#> Loading required package: Biostrings
#> Loading required package: BiocGenerics
#>
#> Attaching package: 'BiocGenerics'
#> The following objects are masked from 'package:stats':
#>
#>     IQR, mad, sd, var, xtabs
#> The following objects are masked from 'package:base':
#>
#>     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
#>     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
#>     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
#>     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
#>     Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
#>     table, tapply, union, unique, unsplit, which.max, which.min
#> Loading required package: S4Vectors
#> Loading required package: stats4
#>
#> Attaching package: 'S4Vectors'
#> The following objects are masked from 'package:base':
#>
#>     expand.grid, I, unname
#> Loading required package: IRanges
#> Loading required package: XVector
#> Loading required package: GenomeInfoDb
#>
#> Attaching package: 'Biostrings'
#> The following object is masked from 'package:ape':
#>
#>     complement
#> The following object is masked from 'package:base':
#>
#>     strsplit
#> Loading required package: DECIPHER
#> Loading required package: RSQLite
#> Loading required package: parallel
#> Loading required package: reshape2
#> Loading required package: phangorn
#> Loading required package: sangerseqR
#> Loading required package: gridExtra
#>
#> Attaching package: 'gridExtra'
#> The following object is masked from 'package:BiocGenerics':
#>
#>     combine
#> Loading required package: shiny
#> Loading required package: shinydashboard
#>
#> Attaching package: 'shinydashboard'
#> The following object is masked from 'package:graphics':
#>
#>     box
#> Loading required package: shinyjs
#>
#> Attaching package: 'shinyjs'
#> The following object is masked from 'package:shiny':
#>
#>     runExample
#> The following object is masked from 'package:RSQLite':
#>
#>     show
#> The following object is masked from 'package:Biostrings':
#>
#>     show
#> The following object is masked from 'package:GenomeInfoDb':
#>
#>     show
#> The following object is masked from 'package:XVector':
#>
#>     show
#> The following object is masked from 'package:IRanges':
#>
#>     show
#> The following object is masked from 'package:S4Vectors':
#>
#>     show
#> The following object is masked from 'package:stats4':
#>
#>     show
#> The following objects are masked from 'package:methods':
#>
#>     removeClass, show
#> Loading required package: data.table
#>
#> Attaching package: 'data.table'
#> The following objects are masked from 'package:reshape2':
#>
#>     dcast, melt
#> The following object is masked from 'package:IRanges':
#>
#>     shift
#> The following objects are masked from 'package:S4Vectors':
#>
#>     first, second
#> Loading required package: plotly
#> Loading required package: ggplot2
#>
#> Attaching package: 'plotly'
#> The following object is masked from 'package:ggplot2':
#>
#>     last_plot
#> The following object is masked from 'package:phangorn':
#>
#>     add_boxplot
#> The following object is masked from 'package:XVector':
#>
#>     slice
#> The following object is masked from 'package:IRanges':
#>
#>     slice
#> The following object is masked from 'package:S4Vectors':
#>
#>     rename
#> The following object is masked from 'package:stats':
#>
#>     filter
#> The following object is masked from 'package:graphics':
#>
#>     layout
#> Loading required package: DT
#>
#> Attaching package: 'DT'
#> The following objects are masked from 'package:shiny':
#>
#>     dataTableOutput, renderDataTable
#> Loading required package: zeallot
#> Loading required package: excelR
#> Loading required package: shinycssloaders
#> Loading required package: ggdendro
#> Loading required package: shinyWidgets
#>
#> Attaching package: 'shinyWidgets'
#> The following object is masked from 'package:shinyjs':
#>
#>     alert
#> Loading required package: openxlsx
#> Loading required package: tools
#> Loading required package: rmarkdown
#> Loading required package: knitr
#> Loading required package: seqinr
#>
#> Attaching package: 'seqinr'
#> The following object is masked from 'package:shiny':
#>
#>     a
#> The following object is masked from 'package:sangerseqR':
#>
#>     read.abif
#> The following object is masked from 'package:Biostrings':
#>
#>     translate
#> The following objects are masked from 'package:ape':
#>
#>     as.alignment, consensus
#> Loading required package: BiocStyle
#>
#> Attaching package: 'BiocStyle'
#> The following objects are masked from 'package:rmarkdown':
#>
#>     html_document, md_document, pdf_document
#> The following object is masked from 'package:shiny':
#>
#>     markdown
#> Loading required package: logger
#> Welcome to sangeranalyseR
data("sangerAlignmentData")
data("sangerContigData")
data("sangerReadFData")

rawDataDir <- system.file("extdata", package = "sangeranalyseR")
fastaFN <- file.path(rawDataDir, "fasta", "SangerAlignment", "Sanger_all_reads.fa")
my_sangerAlignmentFa <- new("SangerAlignment",
                            inputSource          = "FASTA",
                            processMethod        = "REGEX",
                            FASTA_File           = fastaFN,
                            REGEX_SuffixForward  = "_[0-9]*_F$",
                        REGEX_SuffixReverse  = "_[0-9]*_R$",
                            refAminoAcidSeq      = "SRQWLFSTNHKDIGTLYFIFGAWAGMVGTSLSILIRAELGHPGALIGDDQIYNVIVTAHAFIMIFFMVMPIMIGGFGNWLVPLMLGAPDMAFPRMNNMSFWLLPPALSLLLVSSMVENGAGTGWTVYPPLSAGIAHGGASVDLAIFSLHLAGISSILGAVNFITTVINMRSTGISLDRMPLFVWSVVITALLLLLSLPVLAGAITMLLTDRNLNTSFFDPAGGGDPILYQHLFWFFGHPEVYILILPGFGMISHIISQESGKKETFGSLGMIYAMLAIGLLGFIVWAHHMFTVGMDVDTRAYFTSATMIIAVPTGIKIFSWLATLHGTQLSYSPAILWALGFVFLFTVGGLTGVVLANSSVDIILHDTYYVVAHFHYVLSMGAVFAIMAGFIHWYPLFTGLTLNNKWLKSHFIIMFIGVNLTFFPQHFLGLAGMPRRYSDYPDAYTTWNIVSTIGSTISLLGILFFFFIIWESLVSQRQVIYPIQLNSSIEWYQNTPPAEHSYSELPLLTN",
                            minReadsNum           = 2,
                            minReadLength         = 20,
                            minFractionCall       = 0.5,
                            maxFractionLost       = 0.5,
                            geneticCode           = GENETIC_CODE,
                            acceptStopCodons      = TRUE,
                            readingFrame          = 1,
                            processorsNum         = 2)
#> Assessing frameshifts in nucleotide sequences:
#> ================================================================================
#>
#> Time difference of 0.05 secs
#> Assessing frameshifts in nucleotide sequences:
#> ================================================================================
#>
#> Time difference of 0.03 secs
#> Assessing frameshifts in nucleotide sequences:
#> ================================================================================
#>
#> Time difference of 0.03 secs
#> Assessing frameshifts in nucleotide sequences:
#> ================================================================================
#>
#> Time difference of 0.03 secs

writeFasta(my_sangerAlignmentFa,
           outputDir         = tempdir(),
           compress          = FALSE,
           compression_level = NA,
           selection         = "all")

generateReport(my_sangerAlignmentFa,
               outputDir           = tempdir(),
               includeSangerRead   = TRUE,
               includeSangerContig = TRUE)
#> processing file: SangerRead_Report_fasta.Rmd
#> output file: SangerRead_Report_fasta.knit.md
#> /usr/local/bin/pandoc +RTS -K512m -RTS SangerRead_Report_fasta.knit.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /private/var/folders/gf/3p_ynkts411bs238rtw3y0b40000gn/T/Rtmplj2Uy8/SangerAlignment/Achl_ACHLO006-09/Sanger_all_reads/SangerRead_Report_fasta.html --lua-filter /Library/Frameworks/R.framework/Versions/4.2/Resources/library/rmarkdown/rmarkdown/lua/pagebreak.lua --lua-filter /Library/Frameworks/R.framework/Versions/4.2/Resources/library/rmarkdown/rmarkdown/lua/latex-div.lua --embed-resources --standalone --variable bs3=TRUE --section-divs --table-of-contents --toc-depth 3 --variable toc_float=1 --variable toc_selectors=h1,h2,h3 --variable toc_collapsed=1 --variable toc_smooth_scroll=1 --variable toc_print=1 --template /Library/Frameworks/R.framework/Versions/4.2/Resources/library/rmarkdown/rmd/h/default.html --no-highlight --variable highlightjs=1 --variable theme=bootstrap --mathjax --variable 'mathjax-url=https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML' --include-in-header /var/folders/gf/3p_ynkts411bs238rtw3y0b40000gn/T//Rtmplj2Uy8/rmarkdown-str3926187706b9.html
#>
#> Output created: /private/var/folders/gf/3p_ynkts411bs238rtw3y0b40000gn/T/RtmplI8nQX/SangerAlignment/SangerAlignment_Report.html
#> [1] "/private/var/folders/gf/3p_ynkts411bs238rtw3y0b40000gn/T/RtmplI8nQX/SangerAlignment/SangerAlignment_Report.html"

Created on 2023-03-02 with reprex v2.0.2

Example from the docs:
library(sangeranalyseR)
#> Loading required package: stringr
#> Loading required package: ape
#> Loading required package: Biostrings
#> Loading required package: BiocGenerics
#>
#> Attaching package: 'BiocGenerics'
#> The following objects are masked from 'package:stats':
#>
#>     IQR, mad, sd, var, xtabs
#> The following objects are masked from 'package:base':
#>
#>     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
#>     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
#>     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
#>     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
#>     Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
#>     table, tapply, union, unique, unsplit, which.max, which.min
#> Loading required package: S4Vectors
#> Loading required package: stats4
#>
#> Attaching package: 'S4Vectors'
#> The following objects are masked from 'package:base':
#>
#>     expand.grid, I, unname
#> Loading required package: IRanges
#> Loading required package: XVector
#> Loading required package: GenomeInfoDb
#>
#> Attaching package: 'Biostrings'
#> The following object is masked from 'package:ape':
#>
#>     complement
#> The following object is masked from 'package:base':
#>
#>     strsplit
#> Loading required package: DECIPHER
#> Loading required package: RSQLite
#> Loading required package: parallel
#> Loading required package: reshape2
#> Loading required package: phangorn
#> Loading required package: sangerseqR
#> Loading required package: gridExtra
#>
#> Attaching package: 'gridExtra'
#> The following object is masked from 'package:BiocGenerics':
#>
#>     combine
#> Loading required package: shiny
#> Loading required package: shinydashboard
#>
#> Attaching package: 'shinydashboard'
#> The following object is masked from 'package:graphics':
#>
#>     box
#> Loading required package: shinyjs
#>
#> Attaching package: 'shinyjs'
#> The following object is masked from 'package:shiny':
#>
#>     runExample
#> The following object is masked from 'package:RSQLite':
#>
#>     show
#> The following object is masked from 'package:Biostrings':
#>
#>     show
#> The following object is masked from 'package:GenomeInfoDb':
#>
#>     show
#> The following object is masked from 'package:XVector':
#>
#>     show
#> The following object is masked from 'package:IRanges':
#>
#>     show
#> The following object is masked from 'package:S4Vectors':
#>
#>     show
#> The following object is masked from 'package:stats4':
#>
#>     show
#> The following objects are masked from 'package:methods':
#>
#>     removeClass, show
#> Loading required package: data.table
#>
#> Attaching package: 'data.table'
#> The following objects are masked from 'package:reshape2':
#>
#>     dcast, melt
#> The following object is masked from 'package:IRanges':
#>
#>     shift
#> The following objects are masked from 'package:S4Vectors':
#>
#>     first, second
#> Loading required package: plotly
#> Loading required package: ggplot2
#>
#> Attaching package: 'plotly'
#> The following object is masked from 'package:ggplot2':
#>
#>     last_plot
#> The following object is masked from 'package:phangorn':
#>
#>     add_boxplot
#> The following object is masked from 'package:XVector':
#>
#>     slice
#> The following object is masked from 'package:IRanges':
#>
#>     slice
#> The following object is masked from 'package:S4Vectors':
#>
#>     rename
#> The following object is masked from 'package:stats':
#>
#>     filter
#> The following object is masked from 'package:graphics':
#>
#>     layout
#> Loading required package: DT
#>
#> Attaching package: 'DT'
#> The following objects are masked from 'package:shiny':
#>
#>     dataTableOutput, renderDataTable
#> Loading required package: zeallot
#> Loading required package: excelR
#> Loading required package: shinycssloaders
#> Loading required package: ggdendro
#> Loading required package: shinyWidgets
#>
#> Attaching package: 'shinyWidgets'
#> The following object is masked from 'package:shinyjs':
#>
#>     alert
#> Loading required package: openxlsx
#> Loading required package: tools
#> Loading required package: rmarkdown
#> Loading required package: knitr
#> Loading required package: seqinr
#>
#> Attaching package: 'seqinr'
#> The following object is masked from 'package:shiny':
#>
#>     a
#> The following object is masked from 'package:sangerseqR':
#>
#>     read.abif
#> The following object is masked from 'package:Biostrings':
#>
#>     translate
#> The following objects are masked from 'package:ape':
#>
#>     as.alignment, consensus
#> Loading required package: BiocStyle
#>
#> Attaching package: 'BiocStyle'
#> The following objects are masked from 'package:rmarkdown':
#>
#>     html_document, md_document, pdf_document
#> The following object is masked from 'package:shiny':
#>
#>     markdown
#> Loading required package: logger
#> Welcome to sangeranalyseR
rawDataDir <- system.file("extdata", package = "sangeranalyseR")
parentDir <- file.path(rawDataDir, 'Allolobophora_chlorotica', 'ACHLO')

ACHLO_contigs <- SangerAlignment(ABIF_Directory     = parentDir,
                                 REGEX_SuffixForward = "_[0-9]*_F.ab1$",
                             REGEX_SuffixReverse = "_[0-9]*_R.ab1$")

launchApp(ACHLO_contigs)
writeFasta(ACHLO_contigs)
generateReport(ACHLO_contigs)

Created on 2023-03-02 with reprex v2.0.2

My sessionInfo output:
sessionInfo()
#> R version 4.2.2 (2022-10-31)
#> Platform: x86_64-apple-darwin17.0 (64-bit)
#> Running under: macOS Big Sur ... 10.16
#>
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
#>
#> locale:
#> [1] en_AU.UTF-8/en_AU.UTF-8/en_AU.UTF-8/C/en_AU.UTF-8/en_AU.UTF-8
#>
#> attached base packages:
#>  [1] tools     parallel  stats4    stats     graphics  grDevices utils
#>  [8] datasets  methods   base
#>
#> other attached packages:
#>  [1] sangeranalyseR_1.8.0  logger_0.2.2          BiocStyle_2.26.0
#>  [4] seqinr_4.2-23         knitr_1.42            rmarkdown_2.20
#>  [7] openxlsx_4.2.5.2      shinyWidgets_0.7.6    ggdendro_0.1.23
#> [10] shinycssloaders_1.0.0 excelR_0.4.0          zeallot_0.1.0
#> [13] DT_0.27               plotly_4.10.1         ggplot2_3.4.1
#> [16] data.table_1.14.8     shinyjs_2.1.0         shinydashboard_0.7.2
#> [19] shiny_1.7.4           gridExtra_2.3         sangerseqR_1.34.0
#> [22] phangorn_2.11.1       reshape2_1.4.4        DECIPHER_2.26.0
#> [25] RSQLite_2.3.0         Biostrings_2.66.0     GenomeInfoDb_1.34.9
#> [28] XVector_0.38.0        IRanges_2.32.0        S4Vectors_0.36.2
#> [31] BiocGenerics_0.44.0   ape_5.7               stringr_1.5.0
#>
#> loaded via a namespace (and not attached):
#>  [1] nlme_3.1-162           bitops_1.0-7           fs_1.6.1
#>  [4] bit64_4.0.5            httr_1.4.5             R.cache_0.16.0
#>  [7] bslib_0.4.2            utf8_1.2.3             R6_2.5.1
#> [10] DBI_1.1.3              lazyeval_0.2.2         colorspace_2.1-0
#> [13] ade4_1.7-22            withr_2.5.0            tidyselect_1.2.0
#> [16] bit_4.0.5              compiler_4.2.2         cli_3.6.0
#> [19] sass_0.4.5             scales_1.2.1           quadprog_1.5-8
#> [22] digest_0.6.31          R.utils_2.12.2         pkgconfig_2.0.3
#> [25] htmltools_0.5.4        styler_1.9.0           fastmap_1.1.1
#> [28] htmlwidgets_1.6.1      rlang_1.0.6            rstudioapi_0.14
#> [31] jquerylib_0.1.4        generics_0.1.3         jsonlite_1.8.4
#> [34] zip_2.2.2              dplyr_1.1.0            R.oo_1.25.0
#> [37] RCurl_1.98-1.10        magrittr_2.0.3         GenomeInfoDbData_1.2.9
#> [40] Matrix_1.5-3           Rcpp_1.0.10            munsell_0.5.0
#> [43] fansi_1.0.4            lifecycle_1.0.3        R.methodsS3_1.8.2
#> [46] stringi_1.7.12         yaml_2.3.7             MASS_7.3-58.2
#> [49] zlibbioc_1.44.0        plyr_1.8.8             grid_4.2.2
#> [52] blob_1.2.3             promises_1.2.0.1       crayon_1.5.2
#> [55] lattice_0.20-45        pillar_1.8.1           igraph_1.4.1
#> [58] codetools_0.2-19       fastmatch_1.1-3        reprex_2.0.2
#> [61] glue_1.6.2             evaluate_0.20          BiocManager_1.30.20
#> [64] vctrs_0.5.2            httpuv_1.6.9           gtable_0.3.1
#> [67] purrr_1.0.1            tidyr_1.3.0            cachem_1.0.7
#> [70] xfun_0.37              mime_0.12              xtable_1.8-4
#> [73] later_1.3.0            viridisLite_0.4.1      tibble_3.1.8
#> [76] memoise_2.0.1          ellipsis_0.3.2

I used your code and ran the launchapp() function on the contigs. Below is the error code. Listening onhttp://127.0.0.1:6825
Warning: Error in value[[3L]]: Couldn't normalize path in addResourcePath, with arguments: prefix = 'AdminLTE-2.0.6'; directoryPath = 'F:/biocbuild/bbs-3.16-bioc/R/library/shinydashboard/AdminLTE'
3: runApp
2: print.shiny.appobj
1:

Hi, Thank you for responding back. I am using Windows 10 R version 4.2.2. I edited my question to include the sessioninfo() for the package. I also recently commented the error I encountered from using your code. It's probably because my version of R is not updated because its Windows.

Another comment: I tried updating the Bioconductor package but I am using the newest version of it. However, when I used 'BiocManager::install(version='devel')', I ran across this error: Error: Bioconductor version '3.17' requires R version '4.3'; use
BiocManager::install(version = '3.16') with R version 4.2; see
https://bioconductor.org/install'. Is it because I am using Windows? If so, how could I use the example code in another program?

Another comment: I knitted your code again and this was the error I encountered:
`Error in parse_block(g[-1], g[1], params.src, markdown_mode) :
Duplicate chunk label 'unnamed-chunk-1', which has been used for the chunk:
library(reprex)
library(sangeranalyseR)
data("sangerAlignmentData")
data("sangerContigData")
data("sangerReadFData")
Calls:  ... process_file -> split_file -> lapply -> FUN -> parse_block
Execution halted'

Looks like you my have a few issues to deal with. One potential issue is your shiny package install. First thing I'd do is, in a new session, reinstall shiny (install.packages("shiny")) and run the example code from https://shiny.rstudio.com/gallery/plot-interaction-basic.html, then any issues with shiny and related packages should be able to be identified and solved.


text:
LAMP primers design software
What software do you use to design LAMP reaction primers for alignment?  So far I have been using Primer Explorer v5 but I am looking for some alternatives.
Alignment includes some species from the same virus family. So I am looking for a set of primers that will capture species from the given alignment.

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.

Welcome to the site. What is the depth of study and degeneracy, e.g. within a species, genera, or family ... or beyond? This is kinda important to identify which algorithm.

Thanks, woooaa thats difficult. I might have misunderstood this but if its an RNA virus diagnosis at the familial level that is not trivial. The mutation rate of RNA viruses is very high - at familial level, it's impossibly high excepting for the use of nested-PCR strategies.


text:
Non-ACGTN alternate alleles in VCF
I was trying to left normalize my vcf file using following command:
bcftools norm -f reference.fasta -m -any input.vcf > output.norm.vcf
However this gives the following error:
Non-ACGTN alternate allele at 1:278165 .. VCF:'-3AGG'
It looks like the variant caller returned all deletions and insertions like this:
1   5387285 .   C   +2TG
1   5387316 .   G   -8TGTGTATC

This file format was produced by SCCaller 2.0.0, there are several unconventional output fields produced by this caller but but the repository isn't maintained.
Is this an issue of the vcf format or does this means that there is a non-ACGTN character on position 1:278165 in the fasta reference file? I don't see anyting out of the ordinary at this position using IGV.
What would be the best way forward?
I tried the following command to extract non-ACGTN allels but it outputs all alleles:
bcftools view -i 'ALT !~ "^[ACGTN]+$"' input.vcf.gz > errors.vcf
EDIT: including several lines of vcf file
CHROM   POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  SAMPLE1
1   89946   .   A   T   68  .   NS=1    GT:SO:AD:BI:GQ:PL   0/1:False:3,18:0.973:68:273,164,49,117
1   286158  .   A   +2TG    78  .   NS=1    GT:SO:AD:BI:GQ:PL   0/0:NA:11,10:0.982:78:154,97,175,456
1   900119  .   A   G   150 .   NS=1    GT:SO:AD:BI:GQ:PL   0/0:False:17,14:0.995:150:214,136,321,657
1   932613  .   G   -4TTTC  62  .   NS=1    GT:SO:AD:BI:GQ:PL   0/1:NA:4,18:0.998:62:274,165,103,173
1   1168310 .   T   +2CA    35  .   NS=1    GT:SO:AD:BI:GQ:PL   0/0:NA:18,5:0.985:35:77,55,91,712```

That isn't a valid ALT field according to the VCF standard:

ALT — alternate base(s): Comma-separated list of alternate non-reference alleles. These alleles do not have
to be called in any of the samples. Options are base Strings made up of the bases A,C,G,T,N (case insensitive)
or the ‘*’ symbol (allele missing due to overlapping deletion) or a MISSING value ‘.’ (no variant) or an
angle-bracketed ID String (“”) or a breakend replacement string as described in Section 5.4. If there
are no alternative alleles, then the MISSING value must be used. In other words, the ALT field must be a
symbolic allele, or a breakend replacement string, or match the regular expression ^([ACGTNacgtn]+|*|.)$.
Tools processing VCF files are not required to preserve case in the allele String, except for IDs, which are case
sensitive. (String; no whitespace, commas, or angle-brackets are permitted in the ID String itself)

So yes, you can have things other than ACTG, N most commonly or . or *. And you can also have things like <DEL> and <DUP> for structural variants. But a + or a digit make no sense and are not allowed so this is not a valid VCF file.
So, filter them out. You can use something like this:
zcat input.vcf.gz | awk -F'\t' '/^#/ || $5 ~ /^[ACGTNacgtn]+|\*|\.)$/' |
    bgzip > fixed.vcf.gz

I would also report this as a bug to the developers of your variant caller.

Thank you for your comment. I will remove them for now, but is there a safe way to reformat those variants as insertions or deletions in the correct format?

@Macintosh I have absolutely no clue what those things are supposed to be. If I were you, I would use a proper variant caller that knows how to produce valid VCF files. What did you use?

SCCaller 2.0.0.

@Macintosh well that looks like an abandoned project that was last touched 4 years ago (https://github.com/biosinodx/SCcaller). If you absolutely must use that, you need to go through the tool's documentation, understand what kind of variant represenation format they are using and then you could post a new question askiung how to convert those. I mean, I think that C   +2TG should be C  TGTG but it might also need to be C CTGTG depending on how they work. And G   -8TGTGTATC should be... too long to write here, but 8 x TGTGTATC in the REF field and G in the ALT.

@terdon, in the examples given in this question, the number in the ALT is the length of the string of bases quoted, not the multiplicity, i.e. -8TGTGTATC represents the 8 bases TGTGTATC being deleted, and +2TG would be the 2 bases TG being added.
It's still not standard VCF, but this might give a clue as to how to fix the REF/ALT to convert it to the correct format.

Huh, thanks @user36196, that's completely new to me. Is this a known format then? Is it mentioned in the tool's docs? Is there some benefit in this? If I want to know the length of a string, I can measure it, it seems odd to include it.

I would personally use plink2 for this. Nice and simple and unerstandable by all.
plink2 --vcf in.vcf --snps-only just-acgt --export vcf --out out.vcf


text:
eggNOG API, how's it work?
Background eggNOG uses seed orthologues which correspond to its own internal database, so for example J421_1876 can be placed in its online search engine here. This will recover the function and species classification. Brief output

Carbohydrate transport and metabolism gluconolactonase activity

.. and lots of associated species. Thats useful information but not very amenable to coded solutions.
Aim I want to replicate this via its API, i.e. to be placed in a pipeline.
Thus eggNOG API is here - but I can't figure out, here's the instructions

http://eggnogapi5.embl.de/nog_data/[dataformat]/[attributes]/[nogname]

Question Does anyone know how to use the API to interrogate using the seed orthologue code?
There are alternatives which involve LOADS of hard drive space (not very attractive).


text:
Filter phyloseq sequence table
I have phyloseq object where sequence table were builded by using the function mergeSequenceTables from DADA2 package, to merge ten sequence tables from different runs.
The sample_data consist of two columns 1)containing sample names same as sequence table column names and 2)study name from which samples comes from.
Using this information I want to subset sequences that have  average relative abundance > 0.1% at least in 3 studies.
I tried this but it doesn't consider information about studies.
prevalenceDF = data.frame(Prevalence = colSums(otu_table(ps.trans) > 0),
                      TotalAbundance = colSums(otu_table(ps.trans)),
                      tax_table(ps.trans)
)

summary_prevalence <- plyr::ddply(prevalenceDF, "Genus", function(df){ data.frame(mean_prevalence = mean(df$Prevalence),
             totalAbundance = sum(df$TotalAbundance, na.rm = T), stringsAsFactors = F) })

keepPhyla <- summary_prevalence$Phylum[summary_prevalence$totalAbundance/sum(summary_prevalence$totalAbundance)
>= 0.001]
ps.1 = subset_taxa(ps, Phylum %in% keepPhyla)

Does anyone have any pointers haw can I count how many studies containing each genera and add as column to prevalenceDF data frame or other way to do filtering?


text:
Transcript vs Primary transcript on phytozome
Could someone help me understand what the difference between transcript and primary transcript on phytozme is? For example, this dataset of A.thaliana has "primary transcript CDS" vs CDS.
Off the top of my head primary transcript represents the initial strand of RNA made from DNA and transcripts represent genes post-processing of the initial strand.
The thing that is confusing me is, why is the transcript file larger than the primary transcript file? The transcript file has about 21k more headers than the "primary transcript". The only explanation I can think of is that due to alternate splicing you could end up with a bunch of isoforms that lend more headers to the transcript file. Is that it?

According to a team lead at Phytozome:

Multi-exon genes can produce different transcripts via process of alternative splicing:  inclusion/exclusion of particular exons in the spliced transcript.Designation of one of these transcripts as “primary” will vary based on the genome under consideration:  the primary might be the one with the highest expression under “normal” conditions, it might be the one whose mutations was first studied, or it might simply be the transcript with the longest CDS.

In all cases, Number_of_transcripts >= Number_of_primary_transcripts                    Number_of_primary_transcripts == Number_of_genes


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

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

Possible duplicate of this question


text:
COG Enrichment, the overrepresentation test between annotated transcriptome and a set of annotated gene, eggnog-mapper annotation results
I am dealing with a transcriptome of a species and list of genes selected for a gene enrichment analysis. In order to see what kind of genes in my selected list are overrepresented compared with the transcriptome, I ran eggnog-mapper annotation.
After running eggnog mapper, I extract the list of GO, KEGG and COG terms associated.
From there, I would like to do a Fisher-test to highlight COG categories overrepresented in our querry selected list compared to our reference transcriptome.
My questions:

How do you deal with enrichment, this kind of test for COG categories?
For the eggnog results KEGG, EC ...ETC, for example how to assign each KO number to metabolism, • genetic information processing • environmental information processing • cellular processes • organismal systems • human diseases • Drug development
and draw graphs?

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.

Hi @se-tima, welcome to Bioinformatics Stack Exchange. The Stack Exchange question / answer approach works best for single specific questions, and not so great for general discussion. If you want a more interactive experience, I recommend you try a discussion website like reddit. Otherwise, please try to limit your post to a single question. If you have multiple questions, feel free to create separate posts for each one.


text:
Which non-matching DNA bases clash the least?
I am building primers for LAMP which I know will not match all my target sequences perfectly.
To some extent, having primers with multiple variants at variable sites (degenerate primers) and longer sequences (tagged degenerate primers) will help. However, I can't just increase degeneracy (variability) forever because then I will have no amplification efficiency.
Therefore, especially at the 5' end which is less important from an amplification perspective, I will have some calculated mismatches. So I have this puzzle:
If those positions of interest (closer to the 5' end) are likely to be any base with about the same frequency, which is either a) the base or b) the pair of bases, which are likely to clash the least?
My instinctive reaction is that for a single base, a pyrimidine would fit better, because even if it doesn't match, it's smaller? And for a pair of bases, probably a purine and a pyrimidine? I would appreciate your thoughts!

I'm not sure this is the right place to ask something like this since it is about the wet-lab side of things and not bioinformatics. I'm not closing it since I am not 100% convinced, but I suspect it would be better on [biology.se] instead.

I would just use Primer(3)-Blast here:
https://www.ncbi.nlm.nih.gov/tools/primer-blast/index.cgi?GROUP_TARGET=on
I would avoid the approach you outline because:

Primer design algorithms incorporating degeneracy are very good;
I'm not sure you're right (but could be wrong).

The algorithm behind Primer-Blast is Primer3, and is legendary in the primer design business and of course Blast.
I agree that strong conservation at 3' is key to the design, but based on homology rather than chemistry.
The thermodynamic calculation ... I would suggest that a mismatch with AT is going to less severe than GC. GC is a stronger bond, but I'd let the algorithm sort this out and prioritise using personal judgement:

Conservation by homology at the 3' end.
Conservation by homology for GC regions.

Ultimately, primers are cheap so ordering loads of algorithmic designs is usually the best approach.

Just to add an explanation on the subject of why is this bioinformatics?
One of the issues of primer selection is self-hybridisation within the primer or annealing between the primer pair. This is bad news for a primer - because it might not work.
Algorithm primer design programs are good at spotting and avoiding unwanted hybridisation. When degeneracy is involved this becomes complicated.
What Primer-Blast does is assess all combinations of degeneracy against issues such as self-hybridisation and primer-annealing. Humanely it is not possible to assess all those combinations and that is why an algorithm is needed.

Thank you! As helpful as Primer3 is, it is not exactly optimized for LAMP, which is what I'm toying with right now. The alternative, PrimerExplorer5 was insufficient on its own. But it's good advice.

Many thanks @Laura and thanks for the update on PrimerExplore5, I didn't know about it. A good question all round.


text:
In cellranger, can I change the minimum number of reads for autodetecting chemistry?
I ran cellranger multi on a TotalSeq B dataset and got this error
[error] Pipestance failed. Error log at:
sampleUnderscore2/SC_MULTI_CS/SC_MULTI_CORE/MULTI_CHEMISTRY_DETECTOR/_GEM_WELL_CHEMISTRY_DETECTOR/DETECT_COUNT_CHEMISTRY/fork0/chnk0-u0e6a05f224/_errors

Log message:
There were not enough reads to auto detect the chemistry: Sample CITE_2 in "/scratch/palmer/E222532_PR_Run01_051222/outs/fastq_path/AAAVFT3HV"
Note that you can avoid auto detection by specifying the specific chemistry type and version.
- Minimum number of required reads = 10000
- Number of reads available = 1408

I'm not sure what to specify the chemistry as. Does anyone know what it might be for CITEseq as part of TotalSeq B? Alternatively, can I reduce the minimum number of required reads?
Thank you

The read count looks too low for single cell data. Where are you getting your reads from?

I don't think you can reduce the number of reads required for autodetection, and it would be a bad idea anyway, because the guess could go wrong.
From 10X Q&A:

TotalSeq B is specifically designed to be compatible with the Single Cell 3’ v3 and 3.1 workflows

And from CellRanger support:
--chemistry     Optional. Assay configuration. NOTE: by default the assay configuration is detected automatically, which is the recommended mode. You should only specify chemistry if there is an error in automatic detection. Select one of:

    auto for auto-detection (default),
    threeprime for Single Cell 3′,
    fiveprime for Single Cell 5′,
    SC3Pv2 for Single Cell 3′ v2,
    SC3Pv3 for Single Cell 3′ v3,
    SC3Pv3LT for Single Cell 3′ v3 LT,
    SC3Pv3HT for Single Cell 3′ v3 HT,
    SC5P-PE for Single Cell 5′ paired-end (both R1 and R2 are used for alignment),
    SC5P-R2 for Single Cell 5′ R2-only (where only R2 is used for alignment).
    SC3Pv1 for Single Cell 3′ v1. NOTE: this mode cannot be auto-detected. It must be set explicitly with this option.
    ARC-v1 for analyzing the GEX portion of multiome data. NOTE: this mode cannot be auto-detected.

That leaves more than one possible chemistry option that could be applicable to your data. Check with the people who prepared the library if you don't know which library kit was used.


text:
How to run python request for list of url's with multiple page numbers?
Hi I am trying to get the cancer ontologies (obo_id and label) from EBI-OLS. Earlier I have used the below code to get the obo_id terms and labels
descendants = requests.get(
    "https://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0004992/descendants"
)

As there are more than 1 page so I created total_pages
total_pages = descendants.json()["page"]["totalElements"] // 1000 + 1

cancer_id = []

for i in range(total_pages):
    ## Check the number of id's present in 'totalElements'
    des = {"page": i, "size": descendants.json()["page"]["totalElements"]}
    ## Add the number of 'totalElements' and 'pages' in the new python request for descendants
    cancer_efo = requests.get(
        "https://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0004992/descendants",
        params=des,
    )
    ## Get the number of obo_id's present in the cancer_efo.json
    obo = [obj["obo_id"] for obj in cancer_efo.json()["_embedded"]["terms"]]
    ## Get the number of cancer labels present in the cancer_efo.json
    label = [obj["label"] for obj in cancer_efo.json()["_embedded"]["terms"]]
    data = zip(obo, label)
    cancer_id.extend(data)

The above code worked pretty well with one link
Now the issue is that I want to get ontologies for 37 descendent links and maximum links have more than 1 page. I want to get one file (cancer_id) containing all obo_id's and labels from these all 37 links.
I made a separate file called list.txt (Right now I am copying only 3 here)
list.txt
https://www.ebi.ac.uk/ols/api/ontologies/ncit/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FNCIT_C4861/descendants
https://www.ebi.ac.uk/ols/api/ontologies/ncit/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FNCIT_C3171/descendants
https://www.ebi.ac.uk/ols/api/ontologies/ncit/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FNCIT_C3457/descendants

I tried the below code but it is not working
import requests

with open ('list.txt', 'r') as fin:
    links = fin.readlines()

def requestIt(links):
    cancer_all = requests.get(links)
    return cancer_all.json()

data = list(map(requestIt, links))

obo = obj["obo_id"] for obj in data

The error is
obo = obj["obo_id"] for obj in data
                        ^
SyntaxError: invalid syntax

Okay, its done. What you are looking for is as follows:
import requests
from pathlib import Path

# Test code
def get_all_keys(d):
    for key, value in d.items():
        yield key
        if isinstance(value, dict):
            yield from get_all_keys(value)

def requestIt(links):
    cancer_all = requests.get(links)
    return cancer_all.json()

if __name__ == '__main__':

    path = Path('/PathToDir/test', 'link.html')
    with open (path, 'r') as fin:
        links = fin.readlines()

    data = list(map(requestIt, links))

    # Test code
    for x in get_all_keys(data[0]['_embedded']['terms'][0]):
        print(x)

    # Answer code
    for i,d in enumerate(data):
        print (f'Link {i + 1} ')
        for obo_id in d['_embedded']['terms']:
            print (obo_id['obo_id'])

Output

Link 1
NCIT:C8263
NCIT:C9163
NCIT:C5630

Link 2
NCIT:C7175
NCIT:C68700
NCIT:C68696
NCIT:C68697
NCIT:C68698
NCIT:C68699
NCIT:C9155
NCIT:C9288
NCIT:C9287
NCIT:C9019
NCIT:C9018
NCIT:C36055
NCIT:C36058
NCIT:C36057
NCIT:C36056
NCIT:C38377
NCIT:C122688
NCIT:C156718
NCIT:C156719
NCIT:C169107

Link 3
NCIT:C7056
NCIT:C131911
NCIT:C165799
NCIT:C68684
NCIT:C4998
NCIT:C5089
NCIT:C5095
NCIT:C8154
NCIT:C8604
NCIT:C127840
NCIT:C8489
NCIT:C8852
NCIT:C69144
NCIT:C115349
NCIT:C115367
NCIT:C165240
NCIT:C8926
NCIT:C160230
NCIT:C160231
NCIT:C182035

Where Link 1 is ...
https://www.ebi.ac.uk/ols/api/ontologies/ncit/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FNCIT_C4861/descendants

... and so on.
If the code is listed # Test code that is what is used to interrogate the data structure. This doesn't contribute to your answer.
However, given the way I've done it, it's important to understand the data structure. Thus output of # Test code is at the base of this post. If you need another part of the data structure you simply switch id_odo to whatever else you are looking for (at the base of this post).
There will be an easy to use OOP method based out the data object. I did't want to spend time investigating the object ... the question is specifically about obo_id so it was easier to go for that.
Test code output
iri
lang
description
synonyms
annotation
Contributing_Source
Display_Name
Legacy Concept Name
Neoplastic_Status
Preferred_Name
Semantic_Type
UMLS_CUI
code
label
ontology_name
ontology_prefix
ontology_iri
is_obsolete
term_replaced_by
is_defining_ontology
has_children
is_root
short_form
obo_id  # <- its here
in_subset
obo_definition_citation
obo_xref
obo_synonym
is_preferred_root
_links
self
href
parents
href
ancestors
href
hierarchicalParents
href
hierarchicalAncestors
href
jstree
href
graph
href

@Megha
The requestIt(link) function returns a list of dictionary objects.
So your "data" variable is a list of dictionary objects.
Every item on this list (a dict) has a key named "_embedded" whose value is also a dictionary.
This dictionary (the _embedded) has a key named "terms" which is a list of dictionaries.
The obo_ids you need are inside the dictionaries of the "terms" list.
The problem with your code was that you were trying to access the "obo_id" key on the top level dictionary.
So a simple solution would be like this:
with open ('urls.txt', 'r') as fin:
    links = fin.readlines()
resps = [requestIt(link) for link in links] # this one holds all the data
obo_ids = []
for resp in resps:                          # for every response we got
    terms = resp['_embedded']['terms']      # choose the terms list
    for term in terms:                      # for every dict in this list
        obo_ids.append(term['obo_id'])      # extract obo_id and append it to our answer
print(obo_ids)                              # do stuff with obo_ids

You should probably perform some exception handling inside the requestIt() function to handle broken links etc.

Welcome @punxx and this looks cool

@Megha, I have not looked through the code in too much detail, however.
obo = obj["obo_id"] for obj in data

This isn't formal Python code. The correction would simply be to place this in a comprehension ...
obo = [obj["obo_id"] for obj in data]

... unless you want an iterator.
I think there is a deep issue with obj["obo_id"]. If the syntax still throws an exception could you provide a print out of either data please, or else ....
obj for obj in data:
    print (obj)
    print (obj.keys)

The data data structure for the output of get.requests followed by .json() is kinda key here

Thank you for the answer. Got the same error with obo = [obj["obo_id"] for obj in data]. Then I tried obj for obj in data:     print (obj)     print (obj.keys). Got an error `obj for obj in data:
^
SyntaxError: invalid syntax

Is there any way I can run the code (the cancer_id [ ] one)that worked perfectly for one link in a for loop for all the links?

Thanks @Megha. Did you get a print out for print (obj)? Then an error for print (obj.keys)? The idea is to check the output before the error.

Hi got the same error for both print statements?and there was no output, just this error

Hi @Megha over the weekend I will replicate your code and investigate, the object, json, dictionary ... (whatever it is) output. I should get I definitive answer, but I'm afraid I can't do that during this week due to work stuff.

@Megha cancer_all.json() is an object so its behaviour is very different. obj["obo_id"] is not going to work because you'll need to reference the instance. It's possible something like str(obj)['repo_url'] might work. I'll need to check and examine all instances etc ...


text:
How to calculate frequency of a category with respect to the value in another column?
I was trying to calculate the frequency of disease_present (yes) when smoking status is y (yes) for each group (A, B, C, D)
 The formula is = Number of samples with disease_present (yes) / No of samples with smoking status 'y' * 100
df
Sample_code Group   smoking disease_present
1   A   y   yes
2   B   n
3   C   n   yes
4   D   y
5   B   n   yes
6   A   y   yes
7   C   n
8   D   y
9   A   n   yes
10  B   y
11  D   n   yes
12  C   n
13  A   y   yes
14  B   y   yes

I tried this
df.groupby("Group").count()

This gave the same number in every category

Pandas by default turns empty cells into NaN
I tried to replace NaN with empty strings with df1 = df.fillna(""). I still got the same results

I am not able to reproduce the issue.
However I think this is a crosstab problem:
>>> new_df = df.fillna("no")
>>> new_df
# prints:
    Sample_code Group smoking disease_present
0             1     A       y             yes
1             2     B       n              no
2             3     C       n             yes
3             4     D       y              no
4             5     B       n             yes
5             6     A       y             yes
6             7     C       n              no
7             8     D       y              no
8             9     A       n             yes
9            10     B       y              no
10           11     D       n             yes
11           12     C       n              no
12           13     A       y             yes
13           14     B       y             yes

>>> pd.crosstab([new_df.disease_present, new_df.smoking], new_df.Group)
# prints
Group                    A  B  C  D
disease_present smoking
no              n        0  1  2  0
                y        0  1  0  2
yes             n        1  1  1  1
                y        3  1  0  0


text:
Naming functions in R
I have a set of files based on which I expect to load RDS files like R1-Tn5-NHP-55_NA and R1-Tn5-NHP-42_NA but the code is producing just the name of the 55 files not both 55 and 42.
#standardize N set names
rename_KP_set <- function(sample_name_vector)  {
 new_sample_name_vector <- apply(sample_name_vector,
      function(x){
         if (grepl("42",x[3])){
           c(x[1:3],'55',x[5])
         } else if (grepl("55",x[3])) {
           c(x[1:3],'NHP',x[5])
         } else if (grepl("55",x[3])){
           c(x[1:3],'NHP',x[5])
 } else
c(x[1:3],'55',x[5])
 }
}) %>%
     sapply(paste, collapse='-')
     return(new_sample_name_vector)
 }

How I can get separate rds files having 42 and 55 numbers or for that matter any rds with some number attached to them?
I have already posted this question on Stack overflow on February 27th, 2023 but have got no answers so I am reposting that in Stack exchange bioinformatics in search of an answer. Any help on this issue will be considered useful.

Your question has a few wording issues. For example, what do you mean by "...a set of files based on which I expect to load RDS files..."? Provide sample input, expected output and actual output. You'll probably figure out the problem on your own as you prepare these items.

Please [edit] your question and provide us with input files and expected output so we can reproduce the problem and test our solutions. We cannot work on a problem if we cannot reproduce it.

Can you please explain what your function is trying to do? It's not clear to me from just looking at the function.