text:
How to calculate Allele frequency from vcf file
Organism under investigation is  Plasmodium falciparum. How to calculate the allele frequency for each row?

I tried with this code:
    import pandas as pd

# Read CSV file
# df = pd.read_csv("Pf_M76611.pf7.200_samples - Copy.csv")

df= pd.read_csv("/content/Pf_M76611.pf7.200_samples - Copy.csv")

# Define function to calculate allele frequencies
def calc_allele_freq(row):
    # Get total number of samples
    num_samples = len(row) - 9

    # Get reference allele
    ref_allele = row['REF']

    # Get alternative alleles
    alt_alleles = row['ALT'].split(',')

    # Initialize list to store allele frequencies
    allele_freqs = []

    # Loop over alternative alleles
    for alt_allele in alt_alleles:
        # Initialize count for this allele
        count = 0

        # Loop over samples
        for i in range(9, len(row)):
            # Get genotype for this sample
            gt = row[i].split(':')[0]

            # Count this allele if it is present in the genotype
            if alt_allele in gt:
                count += 1

        # Calculate frequency for this allele
        freq = count / num_samples

        # Add frequency to list
        allele_freqs.append(freq)

    # Add frequency for reference allele
    allele_freqs.insert(0, 1 - sum(allele_freqs))

    # Return allele list and frequency list as a tuple
    allele_list = [ref_allele] + alt_alleles
    freq_list = [round(freq, 1) for freq in allele_freqs]
    return (allele_list, freq_list)

# Apply function to each row of the dataframe
allele_freq_df = df.apply(calc_allele_freq, axis=1, result_type='expand')

# Add columns to dataframe
df['Allele List'] = allele_freq_df[0]
df['Proportion List'] = allele_freq_df[1]

# Select columns for output and drop duplicates
# output_df = df[['#CHROM', 'POS', 'Allele List', 'Proportion List']].drop_duplicates()

# # Convert output dataframe to string and print
# print(output_df.to_string(index=False))

Could you please show that table as text format? It's quite difficult to read on my mobile device.

Coding: first and foremost all operations should be performed directly in pandas. It is much more efficient and alot easier. You don't need all the loops, pandas does that. apply is okay. I would use a very different strategy (below).

The actually contents of the csv are not clear.
This is what I think might have happened: there's   multiple VCFs against the same reference I assume these were loaded and a database join operation performed using the reference genome column 'REF' as the index for the join across all input VCFs. I've placed this here because I'm not 100% sure what the input file represents - if it's not what I think it is - thats how I'd get to the allelic frequency.
Why do I think that?

Without multiple strains the allelic frequency cannot be calculated.
The file states '200 samples'
The presence of multiple alleles in the 'ALT' column.

The problems are:

heterozygosity. This assumes falciparum is in its haploid state (it has a diploid stage). If its the diploid stage the allelic frequency cannot be calculated using this input file because the linkage was discarded (which is really important)
The underlying problem is the number of actual strains isn't explicitly stated .
the frequency of the 'ALT' doesn't appear to be declared
even if the a frequency was declared (but not presented) that must be linked to the 'ALT' column, each allele needs a exact frequency, otherwise population genetics will not work on this data.

First conclusion The allelic frequency is performed across each row using the reference genome as the reference via value_counts() divided by the number of strains and presented in its own column. Somehow the relevant columns are missing - will need to see them.
However Again on the information presented I don't think it can be done because if this represents 200 samples the frequency of a given allele is not declared in this data, or simply the number of 'ALT's across all VCFs. Without accessing that number the allelic frequency is not doable because it's entirely dependent on it.
Possible alternative conclusion Its possible the question needs to be redefined not as allelic frequency - this has an exact definition within population genetics, but genetic diversity, i.e. remove the term 'frequency' altogether. What I think is being calculated is genetic diversity of all samples in the file against a reference genome. If that is correct it would be good to change the question accordingly. These are different mathematical entities.
Whatever operation is required is easy to calculate, but its obtaining that information from the input file and the precise nature of the investigation.


text:
Extract specific nucleotide base position for query and subject from text file
I have an input file looks like this:
> CD630_08510_19
Query  121  TGGCTTGGGTTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  121  TGGCTTGGTTTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180

Query  241  GCAATTGTGCTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  241  GCAATTGTACTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300

Query  361  CCATACTTTATAGGTATGTTTATTGGTATGGTTGGATAG  399
            ||||| ||||| ||||||||||||||||||| |||||||
Sbjct  361  CCATATTTTAT-GGTATGTTTATTGGTATGGCTGGATAG  398

The input file is the result file after running the command line blast.
It's a txt file generated after blast alignment and in txt file header format are like this.
I am trying to extract the specific nucleotide base for position 129, 249, 366, 372 and 392 from query and subject aligned sequence. python and perl solution will be fine as well. thank you.
For example: I need my output to have format below like this:
Seq_name        base_position       query_base  subject_base    all_subject_base
CD630_08510_19  129                 G           T               TAT-C
CD630_08510_19  249                 G           A
CD630_08510_19  366                 C           T
CD630_08510_19  372                 A           -
CD630_08510_19  392                 T           C

CD630_08510_20  129                 G           G               GATCC
CD630_08510_20  249                 G           A
CD630_08510_20  366                 C           T
CD630_08510_20  372                 A           C
CD630_08510_20  392                 T           C

all_subject_base: are the bases from position 129, 249, 366, 372 and 392 of subject_base column. I want to keep this in one column thats why all_subject_base.
I tried to extract this with grep but have not been able to. I am new in bioinformatics so any suggestion would be very helpful.

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.

Have you tried looking for ways to get alternative output formats from blast that would be easier to parse?

I agree with @MaximilianPress given the approach used. The thing about Biopython is this already done - so it's easy

FYI, the BLAST text format output is for reading, not parsing with a computer program, and is subject to change without notice. I don't recommend using it, especially if you want to do similar things in the future.

Will you always have just one alignment in your file or can you have muiltiple? Does our solution need to deal with both cases?

And why aren't your numbers in sequence? Why do we go from 180 as the last position of the first line to 241 as the first position of the second? Why isn't that 181? Is that actually what your data look like?

Usually, when you run your BLAST alignment you specify the output format (e.g. as described in "outfmt" in https://www.ncbi.nlm.nih.gov/books/NBK279684/table/appendices.T.options_common_to_all_blast/). If, for some reason, you can't re-run the BLAST command (e.g. if this output file was given to you, or if you can't access your query or reference sequence/s) you can parse this 'readable' output format but it's not straightforward.
Here is a potential solution using awk:
cat positions.txt
129
249
366
372
392

cat input.txt
> CD630_08510_19
Query  121  TGGCTTGGGTTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  121  TGGCTTGGTTTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180

Query  241  GCAATTGTGCTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  241  GCAATTGTACTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300

Query  361  CCATACTTTATAGGTATGTTTATTGGTATGGTTGGATAG  399
            ||||| ||||| ||||||||||||||||||| |||||||
Sbjct  361  CCATATTTTAT-GGTATGTTTATTGGTATGGCTGGATAG  398
> CD630_08510_20
Query  121  TGGCTTGGATTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  121  TGGCTTGGGTTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180

Query  241  GCAATTGTGCTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  241  GCAATTGTACTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300

Query  361  CCATACTTTATAGGTATGTTTATTGGTATGGTTGGATAG  399
            ||||| ||||| ||||||||||||||||||| |||||||
Sbjct  361  CCATATTTTAT-GGTATGTTTATTGGTATGGCTGGATAG  398

awk 'BEGIN {
    OFS = "\t"
}

NR == FNR {
    a[$0]
    next
}

$1 == ">" {
b[++c] = $2
}

$1 != ">" {
for (i in a) {
    if ($1 == "Query" && i >= $2 && i <= $4) {
        d[++e] = b[c] OFS i OFS substr($3, (i - $2) + 1, 1)
    } else if ($1 == "Sbjct" && i >= $2 && i <= $4) {
        f[++g] = substr($3, (i - $2) + 1, 1)
        h[c] = h[c] substr($3, (i - $2) + 1, 1)
        }
    }
}

END {
    printf "%s\n", "Seq_name\tbase_position\tquery_base\tsubject_base\tall_subject_base"
    for (j = 1; j <= e; j++) {
        printf "%s\t%s\t%s\n", d[j], f[j], (j % length(a) == 1 ? h[int(j / length(a)) + 1] : "")
    }
}' positions.txt input.txt
Seq_name    base_position   query_base  subject_base    all_subject_base
CD630_08510_19  129 G   T   TA-CT
CD630_08510_19  249 G   A
CD630_08510_19  372 A   -
CD630_08510_19  392 T   C
CD630_08510_19  366 C   T
CD630_08510_20  129 A   G   GA-CT
CD630_08510_20  249 G   A
CD630_08510_20  372 A   -
CD630_08510_20  392 T   C
CD630_08510_20  366 C   T

If you have more positions:
cat positions.txt
129
249
366
372
373
374
375
376
392

awk 'BEGIN {
    OFS = "\t"
}

NR == FNR {
    a[$0]
    next
}

$1 == ">" {
b[++c] = $2
}

$1 != ">" {
for (i in a) {
    if ($1 == "Query" && i >= $2 && i <= $4) {
        d[++e] = b[c] OFS i OFS substr($3, (i - $2) + 1, 1)
    } else if ($1 == "Sbjct" && i >= $2 && i <= $4) {
        f[++g] = substr($3, (i - $2) + 1, 1)
        h[c] = h[c] substr($3, (i - $2) + 1, 1)
        }
    }
}

END {
    printf "%s\n", "Seq_name\tbase_position\tquery_base\tsubject_base\tall_subject_base"
    for (j = 1; j <= e; j++) {
        printf "%s\t%s\t%s\n", d[j], f[j], (j % length(a) == 1 ? h[int(j / length(a)) + 1] : "")
    }
}' positions.txt input.txt
Seq_name    base_position   query_base  subject_base    all_subject_base
CD630_08510_19  129 G   T   TA-GGTACT
CD630_08510_19  249 G   A
CD630_08510_19  372 A   -
CD630_08510_19  373 G   G
CD630_08510_19  374 G   G
CD630_08510_19  375 T   T
CD630_08510_19  376 A   A
CD630_08510_19  392 T   C
CD630_08510_19  366 C   T
CD630_08510_20  129 A   G   GA-GGTACT
CD630_08510_20  249 G   A
CD630_08510_20  372 A   -
CD630_08510_20  373 G   G
CD630_08510_20  374 G   G
CD630_08510_20  375 T   T
CD630_08510_20  376 A   A
CD630_08510_20  392 T   C
CD630_08510_20  366 C   T

NB this won't work if you only have a single position in positions.txt, but if that is a potential use-case you can modify the script to suit, e.g. replace (j % length(a) == 1 ? h[int(j / length(a)) + 1] : "") with f[j].

Jared_mamrot, Thank you so much.

source: Stack Exchange

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

sort: votes


text:
RNAseq: Why would using a more complete scaffold reduce DEGs?
I have a set of RNAseq data for a couple of relatively under-studied species.
They're fungi, of the order Hypocreales. The assemblies we're working with are de novo - reference seqs don't actually exist yet for most of our species. What the collaborator provided was updated genome assemblies in 60-90 fragments as opposed to the original 250-300, but obviously that's still not full resolution. As for the gene families, not a whole lot is known... in fact our species' genes are still not well described and we're mostly working off of predicted orthology.
During the analysis, a collaborator created new DNA scaffolds for these species which we'd now like to use.
However when I redid our RNAseq using the existing reads with the new scaffolds, one of our species had the number of differentially-expressed genes drop by anywhere from 20-50% per condition.
As the reads themselves haven't changed, we weren't expecting significant changes in differential gene expression, and indeed for our other species changes have been minimal.
Does anybody know what might cause this kind of change?
Thank you!

You mention RNAseq reads, but a DNA scaffold. How are you mapping the reads? Could you please show the command line that was used for read mapping, and give some more detail about how you calculated differential expression?

Fungi, ain't a kingdom most of us will know much about. The only thing I can guess at is that in bacteria there has resurgence in methods for de novo assembly specifically aimed at reducing the total number of contigs in an assembly. The 'new wave' is summarised in this response here, Improving prokaryotic assembly with other contig/scaffold-level data?. This is potentially very exciting, except the OP reported they didn't work ;-(
A reasonable guess is your colleague has deployed these approaches to reduce the number of contigs. HOWEVER, what is appropriate for bacteria, may not be appropriate for fungi, alternatively it may be an improvement. Both bateria and fungi are small genomes, but beyond that there could notable differences.
If fungi are dominated by multigene families these might be being removed in the new wave de novo assemblies (they are designed for bacteria). The collapse of multigene families absolutely occurs in other eukaryotes with small genomes, but they have to be really big on multigene families (>10% of genes in their genome). If that was true you would probably know about it.
Alternative, this could simply be poor gene identification in the first assembly resulting in an artificial jump in the number of genes.
Solution? Either way you should be able to figure which genes are different between the assemblies and precisely what type of genes are differing, e.g. are these large multi-gene families. Personally, I would pipe the predicted genes from both assemblies through eggNOG. This might be useful, but not really what it's designed for, but it will return the KEGG nos.  More explicit methods based around comparing the RNAseq reads against the assembly via e.g. bedtools or these guys here recommend  htseq. This might be a better way to rephrase the question.
Thus I recommend:

Finding out which de novo methods were used between first and second fungi genome assemblies. Do any of those methods match this post, Improving prokaryotic assembly with other contig/scaffold-level data?
Consider a more explicit genomics study to examine RNAseq reads against the genome, or comparison between predicted no. genes between genomes.
Having said all that ... once you've done your DEG you  identify the KEGG no. and then simply compare the lists directly to find out which genes are different between the assemblies. This is an easy database operation to perform, e.g. via a merge function (separate question).

I'd check point 1 and then go for point 3 - because it's something I easily know how to do. Thats just my opinion, a purist would likely explore this via point 2.
Conclusion What I am certain about is that for a species where very little might be known it is better to perform an investigation rather than make an assumption that might prove later to be wrong.

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


text:
sra toolkit fastq-dump comand not terminating
I ran the following command trying to download a fasta file from NIH but the command doesn't seem to terminate and I get no output file in my repository:
fastq-dump SRR10134980

I have installed SRA Toolkit following the instructions on this page.

Do you get the same behavior when you use fasterq-dump? Note that fasterq-dump is the successor to the older fastq-dump tool.

That's an 11GB file and fastq-dump is not a speedster. Without more info the likely answer is that it keeps downloading which might take several hours on  poor connection. You can get fastq download links from sra-explorer.info -- That is more stable than this tool.


text:
Given a FASTA file with a set of genes, how do I determine that set's conservation among a genetic order?
Given a set of genes in a FASTA-formatted file, for example: exFasta.fa, what is the best way to determine how well conserved the genes in that file are among an order of organisms, for example: Orthoptera?
BLAST is an obvious option, but I am not sure that it is the most efficient way considering there are likely thousands of organisms in that order. OrthoFinder also comes across as a valid option, but I am not sure how to implement such a thing considering the required inputs.
Any advice would be greatly appreciated!

Are these annotated, known genes? Or do you also need to detect/identify the potential homologs in the target species?

Orthoptera (grasshoppers, crickets, locust) will be annotated (mostly). What they want to do is skim the draft genomes of Orthoptera. I like @gringers answer, it's a good one. As a project it would work.

It's possible that MOLE-BLAST could help for you:

MOLE-BLAST is an experimental tool that helps taxonomists find closest database neighbors of submitted query sequences. It computes a multiple sequence alignment (MSA) between the query sequences along with their top BLAST database hits, and generates a phylogenetic tree. Query sequences in the tree are denoted with highlighted node labels. If the input sequences come from different genes or loci, MOLE-BLAST can cluster them and compute an MSA and a phylogenetic tree for each locus separately. This tool may eventually be incorporated into our supported suite of BLAST applications based on the feedback we receive.

I haven't actually used MOLE-BLAST, but I have used Smart BLAST, which does a somewhat similar phylogenetic search (although on a restricted diverse taxonomic set) for protein queries.

thank you! I'll try that.


text:
How can I export all and only the publications cited by a specific article?
I want to do a systematic literature review of genetic variant annotation. As part of this I am interested in all the databases available to annotate genetic variants.
I want to export all the references cited by the article with PubMed ID 35943817 (https://pubmed.ncbi.nlm.nih.gov/35943817/), because many of these references are relevant to my review. I can see that they number 93. I want them in a CSV table with their PubMed ID, Title, etc., similar to CSV tables obtained in a search after clicking "Save citations to file", Format "CSV".

Hi @CristianRiccio, your question is fine, but I am not sure this is 'bioinformatics'. The nearest bioinformatics solution is NCBI's Entrez, often used for sequence retrieval, but will retrieve citations from Pubmed - so if you parse the titles in the paper, e.g. Pdfminer you can use those to obtain this data in Entrez, e.g. within Biopython.

You're posting this as a bioinformatics question, which suggests you want code for this. Given that it's only 93 citations, it sounds like you only need to do this once, and they're all in a standard format in the paper. Is there any reason why copying the references from the paper and manually inserting field separators would not be appropriate?

Something like this will work, it did just now for me.  You might want to extract different info from the records though.
This exports to a table (bibliography.tsv) columns containing:

pmid
authors
title
journal
year

max_attempts = 5
def fetch_pubmed_data(pmid):
    Entrez.email = "you@email.com"  # Replace with your email address
    handle = Entrez.efetch(db="pubmed", id=pmid, rettype="xml", retmode="text")
    record = Entrez.read(handle)
    time.sleep(1)
    return record

def extract_cited_articles(pmid):

    record = fetch_pubmed_data(pmid)
    cited_pmids = []

    record_article = record['PubmedArticle']

    if "PubmedData" in record_article[0] and "ReferenceList" in record_article[0]["PubmedData"]:
        for reference in record_article[0]["PubmedData"]["ReferenceList"][0]["Reference"]:
            if "ArticleIdList" in reference and len(reference["ArticleIdList"]) > 0:
                cited_pmids.append(reference["ArticleIdList"][0])

    return cited_pmids

def extract_bibliographic_info(pmid_list, max_attempts):
    bibliographic_info = {}
    bibliographic_info_ls = []

    for pmid in pmid_list:

        authors = title = journal = year = ""
        record = {}

        for attempt in range(max_attempts):
            try:
                record = fetch_pubmed_data(pmid)
                break

            except:
                time.sleep(5)

        if 'PubmedArticle' in record:
            if len(record['PubmedArticle']) >0:
                article = record['PubmedArticle'][0]["MedlineCitation"]["Article"]

                if "AuthorList" in article:
                    authors = ", ".join([f"{author['LastName']} {author['Initials']}" for author in article["AuthorList"] if 'LastName' in author and 'Initials' in author])
                    title = article["ArticleTitle"]
                    journal = article["Journal"]["Title"]
                    year = article["Journal"]["JournalIssue"]["PubDate"]["Year"]

        print (pmid, authors, title, journal, year)
        bibliographic_info_ls.append([pmid, authors, title, journal, year])

    return bibliographic_info_ls

def main():
    pmid = input("Enter the PubMed article ID: ")
    cited_pmids = extract_cited_articles(pmid)
    bibliographic_info_ls = extract_bibliographic_info(cited_pmids, max_attempts)

    df = pd.DataFrame(bibliographic_info_ls, columns=["pmid", "authors", "title", "journal", "year"])
    df.to_csv("bibliography.tsv", sep="\t", index=False)
    return df

if __name__ == "__main__":
    df = main()

Really amazing code. Thank you very much. I just had to change to change cited_pmids.append(reference["ArticleIdList"][0]) to cited_pmids.append(reference["ArticleIdList"][-1]) to get the PMIDs instead of the PMCIDs! Thank you!


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:
What tool(s) can I use to use a .faa file to feed into JSMol?
JSMol, the molecular modelling software, requires PDB file format to work and render the graphics. Currently, I am writing an annotation software - backend in Perl - that takes in a .faa file (fasta amino acids). Is there any way or tool out there that can easily convert to a PDB file? Many of the tools out there do not seem to give what I need.
My workflow so far:
input.fna -> Bakta -> bakta.faa -> covert tool -> proteins.pdb -> JSMol
Essentially, a tool exists called pdb2fasta that takes a pdb file and converts to a fasta file but can’t seem to find a proper tool that does it the other way around.

Or do you just want an unfolded sequence of amino acids (as atoms and bonds)? In any case, you will need something that generates those coordinates, even if just as a linear 'model' of a chain

Would you happen to have any ideas what I could use to do that? @gilleain

@HarishPrabhakar That is an exceptionally difficult problem for a variety of reasons. there is not an obvious way to translate linear AA sequence into three-dimensional structure. You may have heard of alphafold, which uses neural nets to more or less solve the problem. But you have to understand that this is a stochastic and error-prone prediction!


text:
Nanopore Flongle vs MinION
How do the Oxford Nanopore Flongle flow cells perform compared to the standard MinION flow cells, when it comes to reliability and error rates? I found this comparison article from last year, from where the below graph is taken, but they used the old chemistry. How is the situation today with the new R10.4.1 nanopores and Kit 14 chemistry? Does anyone have first-hand experience, or a good reference?
My use case would mainly be metabarcoding, i.e. amplicon sequencing. If anyone has experience of Flongle flow cells for specifically metabarcoding, I would be happy to hear your thoughts.

For us, R10.4.1 Flongle flow cells are running at reduced yield due to lower numbers of working pores (something like 10-20 pores instead of 30-60 pores for R9.4.1 flow cells). This doesn't really affect the stuff we're using Flongles for, because I make sure we've got plenty of headroom built in when designing experiments.
Here are the yield assumptions I use for flow cells, which are basically the lower ranges that we've seen when using them:

Flongle - 200 Mb
MinION - 5 Gb
PromethION - 50 Gb

For Flongle, this has been fine for us to do sequencing and assembly of 12 plasmids on one flow cell using the rapid barcoding kit, and the improved accuracy from the Kit14 + R10.4.1 combination makes reliable assembly easier (especially where homopolymers are involved). I expect it will work similarly well for amplicon sequencing.
Error rates are pretty similar across platforms, because it's the pore that matters, not the package. Here's a comparison plot of mapped accuracy vs predicted accuracy (from the basecaller) that I did recently on a TrackIt 1kb Plus DNA ladder, with bead cleanup using LFB, run on a PromethION flow cell with the V14 ligation kit (sequenced on a P2 Solo), then called using the new bacterial model. The reads represented in this plot have been filtered to only keep (at most) the longest 2000 matches for each target. The run produced 1,128,671 reads (1.95 Gb) over 1.5 hours:

Here's a comparison plot of mapped accuracy vs predicted accuracy (from the basecaller) for the same DNA ladder sample run on a Flongle flow cell (sequenced on a MinION Mk1b), with bead cleanup using SFB, also taking at most 2000 longest sequences from each mapping. The run produced 465,144 reads (191 Mb) over 24 hours:

The predicted accuracy is taken from the basecaller (i.e. fastq quality string), and is calculated by summing up the quality values after converting to a linear representation, then dividing by the number of bases and converting back to a log representation:
sub getMeanQual {
  my ($qual) = @_;
  if(length($qual) == 0){
    return(0);
  }
  my $qBase = 33;
  my @qspl = split("", $qual);
  my $qualTotal = 0;
  my $qualCount = 0;
  foreach my $qual (@qspl) {
    $qualTotal += 10 ** (-(ord($qual) - $qBase) / 10);
    $qualCount++;
  }
  return(int(-10 * log($qualTotal / $qualCount) / log(10)));
}

According to the trained models that LAST generated, there's not much different between the called reads between Flongle and PromethION flow cells. Here's the PromethION base mismatch matrix (q45 SNPs):
# substitution percent identity: 99.9971
#last -Q 1
#last -t4.38507
#last -a 23
#last -A 19
#last -b 2
#last -B 4
#last -S 1
# score matrix (query letters = columns, reference letters = rows):
       A      C      G      T
A      6    -45    -35    -46
C    -45      6    -45    -44
G    -33    -45      6    -46
T    -46    -45    -46      6

And here's the Flongle base mismatch matrix (q41 SNPs):
# substitution percent identity: 99.9922
#last -Q 1
#last -t4.36483
#last -a 20
#last -A 20
#last -b 4
#last -B 5
#last -S 1
# score matrix (query letters = columns, reference letters = rows):
       A      C      G      T
A      6    -45    -36    -46
C    -45      6    -44    -45
G    -32    -45      6    -46
T    -46    -45    -46      6

So I guess you could say that the Flongle has slightly lower accuracy, but I wouldn't call that difference significant for any practical purposes.

Thank you for the update, I have marked your answer as a solution!


text:
Remove files after process is terminated
In my pipeline nextflow, I create a channel with channel.fromSRA and the channel contains lots of heavy files.fastq.gz. Then I have a first process to unzip files and transform them in files.fasta and a second process to count kmers. Is it possible to remove files.gz after first process and remove files.fasta after the second process. Can someone explain me how to do? They are too heavy to keep them in the work directory. Thanks.
In my pipeline I use kmer-ssr to count kmers present on my genome sample.
First I create a reads_ch = Channel.fromSRA
In this channel there are 12 files (. fastq.gz) of 20 Gb each one.
Then the first process unzips them and transforms them in files.fasta.
Now in my work directory there are 12 more files of 80 Gb each one. So the work directory is now 1200 Gb.
The second process counts kmer present in reads and creates 12 files of 100 Gb each one. The work directory will be 2200 Gb.
The work directory is too heavy. I want to know how I can suppress the current file. fastq.gz after first process

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.

Certainly this will be doable. Could you provide the code? Or else sufficient amounts of code to describe the nested loop? I will provide a link below as typical example of a NextFlow question.

So this is yesterdays question here as an example of a standard NextFlow question on the site.

Look into the "temporary" attribute of the inputs that should help you in the right direction.

Unfortunately, this remains an open issue with Nextflow: https://github.com/nextflow-io/nextflow/issues/452
However, writing large amounts of uncompressed data to disk is usually best avoided. If it can't be (e.g. maybe the file needs to be read multiple times), you'd need to merge the two processes so that the writing of the temporary file and the running of the tool is done in a single process. This would then let you remove the temporary input file as soon as it's no longer required.
Another way would be to pipe the uncompressed and FASTA formatted reads directly into your tool. Not all tools support this functionality obviously, but Kmer-SSR can indeed receive input from STDIN, which completely obviates the need to write any temporary files. Note that by default, Kmer-SSR writes also to STDOUT and this lets you write a compressed output file:
process kmer_ssr {

    cpus 4

    input:
    tuple val(sample), path(fastq)

    output:
    tuple val(sample), path("${sample}.tsv.gz")

    """
    zcat -f "${fastq}" \\
        | sed -n '1~4s/^@/>/p; 2~4p' \\
        | kmer-ssr -t ${task.cpus} \\
        | gzip \\
        > "${sample}.tsv.gz"
    """
}

Since the above uses a pipeline, make sure to add the pipefail to your command. One way is to simply apply it to all processes by adding the following to your nextflow.config:
process {

  shell = [ '/bin/bash', '-euo', 'pipefail' ]
}

From the bash man page:

pipefail
If set, the return value of a pipeline is the value of the last
(rightmost) command to exit with a non-zero status, or zero if all
commands in the pipeline exit successfully. This option is disabled by
default.


text:
Using Multi-Dimensional Scaling (MDS) to produce a vector in order to account for patient bias when constructing DGE lists from RNA-seq datasets in R?
I am currently working on my PhD and as part of my thesis, I intend to analyse gene expression within multiple sclerosis (MS) lesions by looking at RNA-seq datasets on Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/)  , which I have selected already. Part of my analysis is using R to transform the count matrices provided on GEO into DGElists (differential gene expression lists), and I have been using both the edgeR handbook, as well as an R script provided by one of the groups whose dataset I am including in my analysis, to construct my own code (as I have had to teach myself R programming from scratch).
[MS lesions are multiple sclerosis lesions. There are multiple subtypes that are included in this group's analysis - Active, Chronic active, inactive, remyelinating (or repairing) lesions and Normal Appearing White Matter. These classifications relate to the physiology of the lesion - active lesions are undergoing active inflammation, chronic active means it goes through periods of inflammation, inactive means inflammation isn't occurring at the moment, and remyelinating relates to a lesion that is undergoing healing.]
My current issue is as follows: Following along with the script provided by the lab group, they appear to have calculated a "PCA Vector", which they use in the creation of their design matrix. As far as I can tell, PCA, or Principal Component Analysis, is meant to account for the multi-dimensionality of a dataset, essentially transforming the data so that it can be plotted on a 2-dimensional set of axes and the variability of the points shown visually. Upon contacting the group, they said that they used PCA to find the biggest outliers as well as identify the largest patient bias, which they then normalised for.
The group has also provided their data which can be opened in R. When I utilise their "PCA vector" in the construction of my design matrix, I get the same values as the group, however, any attempts at producing this PCA vector by performing PCA using R packages does not give me values that are in any way close to those in the group's DGE lists (also provided).
My question therefore is this: Does anybody have an idea as to how one could create a "PCA vector" - a value of variation between patients - which could then be incorporated into a design matrix in order to account for patient variation? I feel I should mention - patient age, sex and MS lesion type are already accounted for as they are also used in the construction of the design matrix. Would anyone have any experience with this?
I ask that you please be gentle/patient with me - I am still very new to the world of bioinformatics (and eager to learn!) and I apologise if I have not explained my issue very well. I will do my very utmost to clear up any confusion, and I would be incredibly grateful to anyone who can help!
For those interested in seeing the journal article which includes links to the GEO repository which hosts their DGE lists and count matrix, as well as their github containing their R script and data, I have included the link below.
https://actaneurocomms.biomedcentral.com/articles/10.1186/s40478-019-0855-7#availability-of-data-and-materials

Welcome @R_Cres_01! "PCA vectors" is kind of unclear. Maybe they mean principal components (PCs)? Do you know which PC they use, and why? Normally PCA is performed on a subset of the (log-normalized) count matrix, using the top N most variable genes. So if you want to recreate their PCA embedding matrix, i.e. a cell x PC matrix where PCs are ordered by variance explained, you need to know how they transformed the data (what normalization? with or without centering and scaling?), how many genes, how they were selected. Also this may not be a good idea, see https://www.biostars.org/p/379656/

Your question title mentions MDS, (multidimensional scaling) but the question itself only mentions PCA (principal component analysis). These are similar concepts, but not the same. Can you confirm you want to do PCA adjustment (more common) rather than MDS adjustment?


text:
How to calculate cell type percentage in every sample
I have a Seurat object (metadata) with the single R samples consisting of cell types and sample types columns. I am trying to make a table that has a sample and percentage of cell types for each sample so that I can run a t-test across the samples that came out of 3 prime vs multi-omics.
I have a seuratobject meta data with singler_labels (cell labels like chondrocytes, astrocytes, macrophage...), sample names (S0030, S0031, S0032.....), cell id/barcodes, and everything a seurat object has. I want to make a table with samples and percent of cell types per sample just to see what sample has what percent types of cells.
How do I do this?
Thank you!

We can't really help you out without looking at your Seurat object. Can you show us what all the metadata in the object looks like?

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, do you have replicates for each condition you want to test? Based on your short description I assume you have but it would be useful to specifically mention this. Be careful when comparing such data because proportions are extremely variable from run to run

If your cell types and sample names are in separate metadata variables attached to the Seurat object, then you can use table to count up the pairings:
table(pbmc.seurat$singler_labels, pbmc.seurat$sample_name)

This will output a matrix containing these counts

The following package performs this type of analysis and can be directly used on a Seurat object:
paper: propeller: testing for differences in cell type proportions in single cell data
github: Speckle
You mention that you would like to do some statistics so keep in mind that this requires replicated measurements to handle the large variability in proportions between samples. Your metadata should contain:
celltype, sample, group_to_test
once the package is installed you can simply run the propeller function and set group to technology (in your post you mention you want to evaluate 3' vs multi-omics):
propel <- propeller(clusters = object@meta.data$celltype,
                    sample = object@meta.data$sample,
                    group = object@meta.data$technology)

This will return the statistics: t-test for two groups or ANOVA for multiple testing.
The following function will return the transformed proportions that can then be used in limma for example to test for more complex designs or simple sample-wise plotting using ggplot. For transform you will offcourse have to decide which transformation is applicable to your data.
getTransformedProps(clusters = object@meta.data$celltype,
                             sample = object@meta.data$sample, transform="asin")

There are a variety of other packages that allow for more complex analysis and a recent benchmark can be found here


text:
How can I make this Biopython program (to correct erroneous barcodes) run faster, and is there any alternative method?
This question has also been asked on Biostars
I am looking forward to getting a valuable suggestion for a bioinformatic problem.
Background:
Currently, I am performing a de novo whole genome assembly. At the stage of barcode correction, I lost nearly half of all the reads due to erroneous barcodes.
Barcodes, in this context, are 18-nucleotide-long DNA sequences that were added to the DNA insert during NGS library preparation. Now, I have a file that contains all the 18-base barcodes used in my projects. These barcodes are important because they are unique for reads that came from a DNA molecule and therefore helpful in generating long-range information.

While library construction, 18-base molecular barcodes were used for the linked-read generation, and the barcodes have a high incidence of "N" in positions between 12 and 18, as you may see below.

There is a very large occurrence of Ns in positions 12 to 18 of the barcodes and for this reason, a large pool of erroneous barcodes were filtered.
Now, I have been asked by the company to do one of two things, to recover the lost reads:

Either replace the barcodes(barcode fastq file from sequencing core) with the correct ones from a pool of barcodes (This is a text file I received from the company that contains the pool of all 18 base barcode sequences)

Since most "N"s are seen in positions 12 to 18 of the 18-base barcode sequences (barcode fastq file from the sequencing core), replace the 11-mers in the barcode file, and replace them with the 11-mers from the correct barcode file (text file).

I have decided to use Python code to execute the first method.
#Define the path where fastq file is stored.
import os
from Bio.SeqIO.QualityIO import FastqGeneralIterator
path=os.path.join(os.path.expanduser("~"), "Desktop", "SVA1_S1_R2_001.fastq")

#create a new file handle to write the file
path_new=os.path.join(os.path.expanduser("~"), "Desktop", "trimmedfastq.fastq")
path_ust_barcode=os.path.join(os.path.expanduser("~"), "Desktop", "UST40MBarcodes.txt")
write_file=open(path_new, "w")

#open the text file containg the correct pool of barcodes

with open(path_ust_barcode, "r") as fh:
   barcodes = fh.readlines()

handle=open(path, "r")
output = []

#===================================================
# Return the Hamming distance between string1 and string2.
def hamming_distance(string1, string2):
    # Start with a distance of zero, and count up
    distance = 0
    # Loop over the indices of the string
    L = len(string1)
    for i in range(L):
        # Add 1 to the distance if these two characters are not equal
        if string1[i] != string2[i]:
            distance += 1
    # Return the final count of differences
    return distance
#===================================================

#set correct barcode to "None"
corr_barcode=None
#set the minimum hamming distance to "None"
min_distance=19
#open the fastq file

for title, sequence, qual in FastqGeneralIterator(handle):

    # if a perfect match exists, just output as is
    if sequence in barcodes:
        output.append("@%s\n%s+\n%s\n" % (title, corr_barcode, qual))
        continue

    for read in barcodes:

        h_distance=hamming_distance(sequence, read)

        # since if there was a perfect match we won't be here,
        # if distance is 1, we'll never do better
        if h_distance == 1:
            corr_barcode = read
            break

        if  h_distance <= min_distance:
            min_distance=h_distance
            corr_barcode=read

    output.append("@%s\n%s+\n%s\n" % (title, corr_barcode, qual))

with open(path_new, "a") as wf:
    for line in output:
        wf.write(line)

handle.close()
print("Done")

Problem
I would like to know if this is the only Pythonic way to tackle this problem. This method is dead slow. My erroneous barcode fastq file is 30+ Gb, and the UST10x barcode file (the text file that contains all the correct barcode pool) is of size 750Mb. I have split this text file into 1000 smaller files, each about 750kb and having ~40,000 lines per file. It's been the third day and the program is still running without outputting anything. I am looking forward to getting valuable suggestions here, with respect to the code, especially how I can make the process faster.
NB:  I've tried replacing
with open(path_ust_barcode, "r") as fh:
   barcodes = fh.readlines()

with
barcodes = list()
with open(path_ust_barcode, "r") as fh:
   for b in fh:
      barcodes.append(b.rstrip())

Still, the code runs extremely slowly.
I am using an IBM server with a 197GM RAM, and an octa-core processor. In my initial code, it was iterating through all the correct barcodes of the text file (UST10x file), for each erroneous barcode in the fastq file, at the end of every iteration, it opened a text file in append mode and wrote the correct version of fastq reads. Since my server has high memory capacity, I rewrote the code in such a way that it stored all the correct barcodes in the memory, and finally wrote it in a text file, by bringing the write function outside of the loop.
Python is not absolutely necessary. I'm not familiar with C or C++, but that being said, I would really love to get this process done. I'm using an IBM server with Ubuntu (Jammy), and 197GB Ram, and the current local storage has only 150GM free space. if you can provide me suggestions as to how to execute it in C or C++, I can do it.

Comments have been moved to chat; please do not continue the discussion here. Before posting a comment below this one, please review the purposes of comments. Comments that do not request clarification or suggest improvements usually belong as an answer, on [meta], or in [chat]. Comments continuing discussion may be removed.

Try profiling (benchmarking) your code to work out where it is slow.
My other suggestions are:

replace your simple Pure Python hamming_distance function with something optimised like RapidFuzz - see https://github.com/maxbachmann/RapidFuzz

cache the corrections you are making - you have tonnes of RAM so a building a dictionary mapping the fragments with N to the corrected fragment ought to be possible. And then for each new sequence, check if you already solved it and if so reuse the same choice.

Thank you for adding a response. I have posted the same question on Biostar, where you can find the code aligned with your logic. A modified version of the code was suggested by @IanSudbery, and you may take a look at it. https://www.biostars.org/p/9561109/#9561603


text:
calculate the backchain RMSD between two pdb files by pymol in Python
i have two protein pdb files and want to calculate the backchain RMSD between them. As far as I know the GUI of Pymol can use align to calculate the RMSD (but I don't know if it is the backchain RMSD).I want to use the python code to use pymol to calculate the backchain RMSD. I know how to align them but don't know how to calculate the RMSD.
import pymol2
with pymol2.PYMOL() as pymol:
    pymol.cmd.load('foo.pdb', 'foo')
    pymol.cmd.load('bar.pdb', 'bar')
    pymol.cmd.align('foo', 'bar')

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.

The backbone can be specified using a selector for the names of these, which are N, CA (Cα) and C (C-prime).
out = pymol.cmd.align('foo and name N+CA+C and polymer', 'bar and name N+CA+C and polymer')
rmsd, n_atoms, n_cyles, n_rmsd_pre, n_atom_pre, score, n_res = out

If you do not wish the proteins to be moved you can use rms_cur, which is just the rmsd.

Thank you for your answer! It works well! however,i want to further learn related knowledge about pymol like this. Can you provide some web/cook book name?


text:
Arlequin - How to assign samples to populations?
I have converted a VCF file to the arlequin file format. It consists of 20 samples that should belong to three populations, along with a few thousands SNPs. However, when I load the arlequin file into the software, it shows I only have 1 population. Any way to fix this?
Thanks.

Can you please show a few lines from your Arlequin file?

Arlequin format is historically extremely difficult to manually manage. The group leader eventually wrote software to allow easy Arlequin format management and conversion. This is achieve by Population Genetics Data (format) Spider, or PGDSpider by Lisher and Excoffier (2018). All information is available here and the manual is here. I would personally avoid attempting to manually rectify population designation.


text:
how to do manual search on NCBI with biopython and add link to publication
I'm trying to reproduce a complex manual search on NCBI using BioPython, but obviously something is not right. My command line is not working, but I'd especially like a way to have an accurate search where I can select only the genomes related to a published paper. Is there an easy way to do that?
I want to limit my results to only those genomes that i) have an NCBI BioProject entry and ii) that BioProject entry contains a Publications field.
For example, let's say I want all the genomes of Microbacterium, the latest version, from RefSeq and not partial (this is what I tried to write but something is not working) and I would like to add one more filter. Take, among the remaining/selected genomes, only the ones related to published papers.
For this strain, we have a related paper:
https://www.ncbi.nlm.nih.gov/bioproject/PRJNA669461/
For this one, we do not:
https://www.ncbi.nlm.nih.gov/bioproject/PRJNA608230/
and I really would like to not have to check 100 genomes manually, but I don't know how to manage it.
I tried:
query = "Microbacterium[Organism] AND latest refseq[filter] AND (all[filter] NOT
partial[filter])"
handle = Entrez.esearch(term=query, db="Assembly", retmax=600)
ids = Entrez.read(handle)["IdList"]

And after this, I'm just using a for loop to get genomes, but my query is not good and it is not handling the problem of the publication.
Thanks for any help!

Edit : as I need to search in the db Assembly unfortunately the answer with the flag author or date is not working

Entrez has the tags [journal], [volume], [page number], [title], [year] and [author] to identify all sequences from a given publication. You will find a full list of tags here from the drop down menu.
Thus,
query = "Microbacterium[Organism] AND Nature[Journal] AND 528[volume] AND Bai[author] AND 364[page number]"

This will identify all Microbacterium genome sequences associated with Bai et al (2015) Nature 528:364- and is here.
For practical purposes I routinely search using the organism and author tags for their sequence output.

The example given (post-edit) was

Kim KK et al., "Microbacterium aquimaris sp. nov., isolated from seawater.", Int J Syst Evol Microbiol, 2008 Jul;58(Pt 7):1616-20

query = "Microbacterium[Organism] AND 58[volume] AND Kim[author] AND 1616[page number] AND 2008[year]"

I removed the [journal] tag because the abbreviation it uses is a bit unusual and replaced it with [year]. Still works here.

From the comments ... if only published sequences were required then using the range function and [publication date] is one way in ...
query = "Microbacterium[Organism] AND 2000/1/1:2020/1/1[publication date]"

The best strategy is to use the [author] because data production tends to be dominated by a small number of groups.

Comments of Comments
Why isn't searching working. The following search works
query = "Microbacterium[Organism] AND (all[filter] NOT partial[filter])"

The term it doesn't like is latest refseq[filter] if that is removed it will work. refseq works, but it doesn't like latest.

I think the OP wants to limit the results to only those genomes whose NCBI BioProject entry includes a Publications field.

@terdon I don't believe there is a tag [publication], there is [publication date], which is better than [year]. However, it's the same answer, the sequence output of most fields are represented only by a small number of authors and only those authors are searched.

@M__ thanks a lot for your help ! your range idea it's pretty clever, it will already be much better than checking by hand!

Thanks @MarineBergot please remember to "accept" if it's useful. You can always switch "accepted" posts if there's a better way posted later.

@M__ by the way, do you have any idea why my current research is not returning anything?

@MarineBergot I've answered in the comments of comments

@M__ thanks for your answer. I arrived in the same point but i need both, latest and refseq. For me latest is working with organism and no partial but refseq is not working even alone. do you have results with Microbacterium[Organism] AND refseq[Filter]? could it be related to version ?

@MarineBergot its the database. The output of the above search using the nucleotide database is here


text:
How to get a file with the number of reads for several fastq.gz files?
I have generated several FASTQ files and I would like to know the amount of reads for each of them.
I am planning to run FastQC on the files which I know would give me the number of reads per sample but it would be a lot easier for me if I could get one text file with the file name and number of reads for all samples together.
When I check for one file I normally do zcat filename_R1.fastq.gz | wc -l and then divide by 4.
Is there a way to do it for all my files at once and get a text file with the file name and reads at the end?
For example I will have 71 R1.fastq.gz files and I would like to output the amount of reads to a text file. I just meant that I would prefer to use one command instead of having to use the command I mentioned above for each file.

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.

Note that fastq files can have more than 1 line of sequence per entry so dividing by 4 is not robust. Make sure to check your file and see if you have cases of multi-line sequences.

One way would be to use seqkit stats, which can be installed using conda:
conda install -c bioconda seqkit

It has a number of useful options, including ways to parallelize the counting, using -j, --threads (default 4), and to skip errors using -e, --skip-err. The latter can be especially useful if you have a larger number of FASTQ files with millions of reads and you'd like to run the command overnight for example. There's also an option to produce a machine-friendly tabular output, using -T, --tabular which would allow you to parse the output without much effort. I'd be tempted to just write the output to a file using the -o, --out-file option:
seqkit stats -To stats.tsv *.fastq.gz

Example results:
$ cat stats.tsv
file    format  type    num_seqs    sum_len min_len avg_len max_len
PHB_Rep1_R1.fastq.gz    FASTQ   RNA 118571  11857100    100 100.0   100
PHB_Rep1_R2.fastq.gz    FASTQ   RNA 118571  11857100    100 100.0   100
PHB_Rep2_R1.fastq.gz    FASTQ   RNA 144826  14482600    100 100.0   100
PHB_Rep2_R2.fastq.gz    FASTQ   RNA 144826  14482600    100 100.0   100
PHB_Rep3_R1.fastq.gz    FASTQ   RNA 129786  12978600    100 100.0   100
PHB_Rep3_R2.fastq.gz    FASTQ   RNA 129786  12978600    100 100.0   100
PUH_Rep1_R1.fastq.gz    FASTQ   RNA 227392  22739200    100 100.0   100
PUH_Rep1_R2.fastq.gz    FASTQ   RNA 227392  22739200    100 100.0   100
PUH_Rep2_R1.fastq.gz    FASTQ   RNA 162373  16237300    100 100.0   100
PUH_Rep2_R2.fastq.gz    FASTQ   RNA 162373  16237300    100 100.0   100
PUH_Rep3_R1.fastq.gz    FASTQ   RNA 185442  18544200    100 100.0   100
PUH_Rep3_R2.fastq.gz    FASTQ   RNA 185442  18544200    100 100.0   100

And then just parse the output to get what you need, if or when you need it. For your specific requirements, you could just use AWK to skip the header line and print the first and fourth columns:
awk 'BEGIN {FS=OFS="\t" } NR>1 { print $1, $4 }' stats.tsv

Results:
PHB_Rep1_R1.fastq.gz    118571
PHB_Rep1_R2.fastq.gz    118571
PHB_Rep2_R1.fastq.gz    144826
PHB_Rep2_R2.fastq.gz    144826
PHB_Rep3_R1.fastq.gz    129786
PHB_Rep3_R2.fastq.gz    129786
PUH_Rep1_R1.fastq.gz    227392
PUH_Rep1_R2.fastq.gz    227392
PUH_Rep2_R1.fastq.gz    162373
PUH_Rep2_R2.fastq.gz    162373
PUH_Rep3_R1.fastq.gz    185442
PUH_Rep3_R2.fastq.gz    185442

You can replace gunzip -c with zcat.
#USAGE: script.sh /runFolder
for i in $1/*.1.fastq.gz
do
   echo "$i" $(gunzip -c $i | echo `wc -l`/4 | bc -l)
done

#OR in the current folder, do this:
for i in ./*.1.fastq.gz
do
   echo "$i" $(gunzip -c $i | echo `wc -l`/4 | bc -l)
done

Backticks are legacy syntax. Subshells can be nested for example using: echo "${i}: $(echo "scale=2; $(zcat "${i}" | wc -l)/4" | bc)"

Also, I think the OP is looking for a summary text file, so you could just have your loop append to an output file (i.e. ... >> counts.txt)

zcat, like cat, works with multiple file names as arguments, up to the argument limit (on Linux, that's usually at least a thousand), so counting multiple lines in multiple files still works, assuming the fastq files are formatted in the expected way with only four lines per record:
zcat filename_R1.fastq.gz filename_R2.fastq.gz ... | wc -l

If there's a consistent file pattern, that can be used with file name glob wildcards:
zcat *.fastq.gz | wc -l

If you have a file containing a list of files, this can be processed as shell arguments as well:
zcat $(cat list_of_files.txt) | wc -l

For an unlimited number of files, something like GNU parallel can be used to feed file names from piped input. It's also a bit safer for file names containing weird characters in them:
cat list_of_files.txt | parallel zcat {} | wc -l

... but [after checking the question again] given that you want counts for individual files, the file name needs to be added as well:
cat list_of_files.txt | parallel echo {} $(zcat {} | wc -l)

The read function provides an alternate way to do the same thing:
cat list_of_files.txt | while read fileName; do echo ${fileName} $(zcat ${fileName} | wc -l); done

Assuming you are using bash, dividing by 4 can be done via the arithmetic substitution operator, $((X)):
cat list_of_files.txt | parallel echo {} $(( $(zcat {} | wc -l) / 4 ))

Or with read:
cat list_of_files.txt | while read fileName; do echo ${fileName} $(( $(zcat ${fileName} | wc -l) / 4 )); done

And, finally, this can be output to a file using >:
(cat list_of_files.txt | while read fileName; do echo ${fileName} $(( $(zcat ${fileName} | wc -l) / 4 )); done) > fastq_record_counts.txt

You'll need to wrap your commands in single quotes to prevent the dollar variables from expanding. You can also avoid the Useless Use of Cat (UUoC) by redirecting stdin or using the :::: argfiles syntax with GNU parallel. To be safe, remember to always quote your variables.

For echo? I don't think it'll make any difference

Sorry, I mean the parallel commands need to be singled quoted: parallel 'echo {} $(( $(zcat {} | wc -l) / 4 ))' :::: list_of_files.txt

If you're looking for a one-liner,
for i in $(ls *.fastq.gz); do count=$(zcat $i | grep "^@" | wc -l); printf "%b\t%b" $i $count; done > counts.tsv

Note that fastq files can include '@' in the first character of the quality string, so this may produce incorrect results.

Also, NEVER do for i in $(ls *), that is very bad practice, will break on even slightly odd file names and is really not needed since a simple for i in  *.fastq.gz; is safer, more robust and shorter to write.


text:
How to do pathway analysis after scanpy for single cell data after DE analysis?
I have done DE analysis using SCANpy on my single cell data and I have compared each cluster versus all the other clusters. One cluster seems particularly interesting so I wanted to do pathway analysis.
I have already tried to do the pathway analysis but I got confused. I initially followed the instructions from here (https://nbisweden.github.io/workshop-scRNAseq/labs/compiled/scanpy/scanpy_05_dge.html) and did the enrichment and then I also wanted to do GSEA so I used gseapy.prerank. However, then I was trying to find a way to visualize the network so I ended up going to the gseapy documentation.
The steps of the enrichment analysis seemed the same but then I got confused with the GSEA portion. The first tutorial I found used gseapy.prerank but the second one uses gp.gsea. I was trying to understand the difference and it seems that since I already have the DE analysis I should be using the gseapy.prerank right bit I am not sure.
I was wondering if doing the analysis trough KEGG would be easier.
Thank you

Gene set enrichment analysis (GSEA) is a commonly used algorithm for characterizing gene expression changes. However, the currently available tools used to perform GSEA have a limited ability to analyze large datasets, which is particularly problematic for the analysis of single-cell data. To overcome this limitation, we developed a GSEA package in Python (GSEApy), which could efficiently analyze large single-cell datasets.

An easy way to do an enrichment analysis in scanpy is using scanpy.queries.enrich(). This relies on gprofiler, which uses KEGG, GO and other ontologies

The original question was GSEApy vs KEGG. The response below is GSEApy.  The question has been edited and is now a  different question concerning options within GSEApy.

Response to original question
If GSEApy as described in the publication provides the output you require - that is much, much simpler. KEGG will be complex and you'd need its ids.
However, surely GSEA will require to perform the full calculation of fold-change from scratch rather than be fed the output of scanpy? KEGG will require the output of scanpy. GSEApy primary function is fold-change rather identifying the precisely which pathways are regulating the phenotype under analysis.
KEGG will provide a more accurate answer - in my humble opinion. GSEApy is a lot easier because it very much an "all-in-one" solution. If the answer isn't satisfactory you can always switch.

Just to note the GSEApy version I'm referring to is the latest publication.
NOTE There was a deleted post from the OP about an option within GSEA. This is a separate question in the SE format of questions and answers. Thats the way the site works.


text:
How to plot average gene expression in scanpy?
I would like to make a UMAP where the cells are colored by the average expression of the bulk signature genes but I am not confident that I did it correctly. I would like to use scanpy for it.
I did the below:
bulk_de_genes_up_list = bulk_de_genes['Gene'].tolist()
# Subset the data based on the list of genes
adata2 = adata[:, adata.var_names.isin(bulk_de_genes_up_list)]
average_expression = adata2.X.mean(axis=1)
adata2.obs['bulk_de_gene_average'] = average_expression
sc.pl.umap(adata2, color='bulk_de_gene_average', cmap='viridis')

I do get a UMAP as an output but I am not sure if it is done correctly. I am mainly worried about average_expression = adata2.X.mean(axis=1)
Is that the correct way of calculating the mean of the gene expression per cell?
Thank you


text:
what does ALT ID DUP refer to when ALT ID for copy number are already specified?
ALT ID CN2, CN3,CN4,CN5 are regions of elevated copy number relative to the reference, Then what does the ID =DUP separately cover? I am new to these analyses. please help
##ALT=<ID=DUP, Description="Region of elevated copy number relative to the reference">
##ALT=<ID=CNO, Description="Copy number allele: 0 copies">
##ALT=<ID=CN2, Description="Copy number allele: 2 copies">
##ALT=<ID=CN3, Description="Copy number allele: 3 copies">
##ALT=<ID=CN4, Description="Copy number allele: 4 copies">
##ALT=<ID=CN5, Description="Copy number allele: 5 copies">

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.

The details will depend on the specific tool used to generate the VCF. Each tool has its own system and there is not much standardization for these fields. That said, the likeliest situation here is that if the CNV caller cannot determine the copy number, but only that there is a duplication, it just puts DUP. And if it can determine the specific copy number, then it uses one of the specific terms.
Alternatively, if the copy number is beyond those numbers with a specific CNN value (say you have 7 copies, or 65 or whatever), then it would presumably also use a DUP.

Thanks for clarifying @terdon, this output is from Canvas 1.39


text:
Where/How can I download all proteins in the UniProt with signal peptide annotated?
Where/How can I download all proteins in the UniProt with signal peptide annotated? I downloaded the uniref50.xml. However, the annotation information is not contained.
In an XML record of a protein, the signal peptide information is contained in <feature type="signal peptide" evidence="2"> </feature>. The uniref.xml does not have this information. uniprot annotates this feature with signalp mainly.

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

Do you have a reason to expect that the signal peptide should be annotated there? Can you point us to more information about how uniprot annotates such features?


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.

For future reference, you can use bcftools norm with the option --check-ref x to exclude such sites.

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

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

DESeqDataSetFromMatrix(se, design ~ designCol)

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

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


text:
How to calculate the variant allelic fraction (VAF) from VEP output?
This question was also asked on Biostars
How can I calculate the variant allelic fraction from the VEP file for multiple human samples?
I used the ensembl-vep and it is human data.

Please show a) your VEP output and b) describe what variant allelic fraction  means. Thank you.

If you just have a VEP output (i.e. VCF) without raw genotypes then you cannot calculate the allelic fraction.

thanks user438383, how can calculate it if i have the vcf file?

@Joman See: https://www.biostars.org/p/147443/


text:
Misidentified cell annotation from azimuth
I'm new to scRNAseq analysis, and am learning how to use Seurat v5 with lung datasets from public data sources (3 GSEs). my workflow is like this:

Load the datasets,
Filter by nFeature and percent mt
Run azimuth (on those separate datasets)
Merge datasets into one Seurat object
Normalize data, find variable features, scale, and run PCA
Integrate layers with harmony
Find neighbors and clusters

By this step, I found that of the cells annotated by azimuth as epithelial (in level 1), some of them at the finest level of detail are misidentified as immune and endothelial (cd8, cd4, macrophages, EC arterial and so on), even though the amount is very small.
Picture attached.
My questions are :

Was my workflow right?

Why are there some non epithelial cells (in the finest level) identified as epithelial cells (in level1) by azimuth annotation?

What should I do for the problem number two; am I safe to just ignore/remove it for further downstream analysis?

Thank you very much.
Epithelial cells (annotated by azimuth) with marker for epithelial, immune, and endothelial:

Cell count and visualization of proportion of each cells:


text:
Can Cytoscape 3 export network.js as UTF-8?
When I set multibyte string(Japanese words) to a node label and export as Full web application, an exported 'network.js'(in 'web_session/data' directory) is encoded as 'ISO 8859-2'(detected by VSCode automatically). And when I reopen 'network.js' as UTF-8, the node label is not encoded correctly.
Can Cytoscape 3 export network.js as UTF-8?


text:
Calculating Z-score from count data
I'm new to this and I can not find information about this:
i have rna seq- I created count data
now I have 2 groups controls and patients.
now I wanted to calculate z scores from this data
and I wanted to know what is the right way to do it:

do I use the z = (x – μ) / σ directly on the count data
2 The μ is calculated per group or from all samples (controls and patients togetherץ

I have RNA seq of different individuals (this is patients rna seq) with no replicates. the sequencing was done with the same method for all of the samples. I want to look at the differences in expression between both groups and because of that I wanted to calculate z score. in this case do I need to calculate the μ per group?
I have ~800 controls and ~1800 patients; but it's not a clinical trial.

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

Please provide more context around this problem; there's not enough information here to give a good answer. For example, how many genes are you looking at - was this a whole-transcriptome study, or a targeted study? How was the cDNA sequencing data generated?
In particular, you shouldn't be using a straight Z-score for RNA/cDNA sequencing data because it's not a normal curve fit in linear space; it is best done using a properly-modeled fit curve. See here for how to process RNA/cDNA sequencing data using DESeq2:
https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#quick-start
Note that DESeq2 requires biological replicates, but you appear to have them in this case: ~800 control replicates, and ~1800 non-control replicates. Given that DESeq2 is able to handle single-cell data when used by Seurat, there's probably a way to work around scale issues if they exist.
However, DESeq2 also expects that the majority of observed genes are not differentially expressed (as is typical for a whole-transcriptome study), which allows it to model the expected expression distribution in order to identify outliers. It won't work for single genes, and is unlikely to work for targeted panels where the majority of genes are differentially expressed. There's a bit of discussion around modeled dispersion here.

Its assumesd the OP is fitting post normalisation. Z-scores are a standard approach to patient cohorts. One issue is the double transformation, but double transformations happen. This is a full scale stage 2 clinical trial, its alot more complicated than described here. It would be against my contract to help, because this is a major big bucks project UNLESS they have a small trial and are scooping the rest of the data from the public domain

I don't think post-normalisation is a reasonable assumption to make. Adi mentioned, "I'm new to this", and that they have count data (which is not normalized).

Maybe but there must statistical expertise on a project of this size

I agree; this is an XY problem, and therefore a direct answer to the question asked is unlikely to fix the actual problem.

Thank you for all the advises I did try to do it with deseq2. I will try to understand it more deeply before doing the analysis. Thank you.

It's not a clinical trial.

@AdiGershon I've updated my answer a little bit. If you've got 800 controls and 1800 patients, then you do have replicates.

Below is the answer for standard normalisation via the z-score. This is RNAseq data and that gets complicated because you'll need an internal standardisation per sample because variation can occur between runs. What we appear to be talking about here is a standardisation over a standardisation. This needs more information about the experimental conditions of the testing because this could result in an erroneous transformation.
Standardisation through a timeseries per patient remains okay however - its difficult to go wrong on that.
Otherwise I'd wait for the input of @RamRS, he may not be around for Christmas.

You wrote:

2 The μ is calculated per group or from all samples (controls and patients togetherץ

Its per group but if its a longitudinal study its per patient (below). As described in point 2 be very clear about you are standardising and why; thats the key to the calculation. You are standardising for an uncontrolled parameter BUT the group (cohort) has something different about it e.g. it's not the drug you are testing but the conditions of the testing. For example, group 1 and group 2 were tested at different times by different people THEN standarisation is useful.
If its just lets standardise stuff when group 1 and group 2 are really identical except for the thing in the test (e.g. a drug vs a placebo); in that scenario standisation needs very careful thought.

To answer the question directly, if you have one patient with longitudinal data then the mean (μ) and standard deviation (σ) will be for the entire longitudinal group for that given patient. x will be for each individual  longitudinal reading. You tend to remove points that are great than 3σ

What you are doing is standardising that given patient so you can now compare values between patients. For example, the time/date of the readings will often differ for different patients and e.g. maybe environmental temperature plays a role ... that sort of thing.

To be clear in the example above you don't take μ for the entire cohort - in a longitudinal study its per patient. If the study is not longitudinal enough e.g. identically performed cohorts - but each cohort differs in some uncontrolled way to the other cohort then the  μ is for a given cohort, but you need to think carefully what you are standardising.

Hi @AdiGershon if you update yoyr post. This isnt a trivial experiment also all the heterogeneiry in sampling and experimental sensitivity in time. Its impossible to simulataneoulsy sample 2000 patients and i doubt the same person did that. At present a study on that scale requires serious attention. I am not entirely sure its doable via this format

Thank you for you comment. I understand it's more complicated then I thought


text:
How do I group sequences by species with a guide tree in IQ-TREE for -asr?
I have an alignment file of DNA sequences with fasta headers in the format
>species-name_protein-name

and a Newick-format (.nwk) tree of the species. The DNA sequences are aligned, because they differ in length.
I am trying to run
iqtree -s sequence.phy -m MFP -msub nuclear -asr -te species.nwk

but run into two problems:

If I just run this, I end up with the protein names not matching the
guide tree.
If I change the protein names to the species names, I
have duplicates.

How do I get IQ-TREE to link the guide tree and the alignment correctly?
The main goal here is to obtain the -asr output.
Any solutions including using a different software than IQ-TREE to find the asr are welcome.
The species.nwk file contains a tree of (in this case anglerfish) with branch lengths and species names:
((Lophius_americanus:6.16014000,Lophius_piscatorius:6.16014000)'14':69.50036000,(Antennarius_maculatus:46.40414000,(Antennarius_striatus:42.06712000,Brachionichthys_hirsutus:42.06712000)'13':4.33702000)'25':29.25636000);

The first couple of sequences from the DNA sequences are
>Lophius-piscatorius_GP-Orthoebolavirus-restonense
AGAGTGGCTCACTCCCTACCCAAGAAGCCTCAGTTGAAGCGGGTACGGCGAGAAGATTCCGATATCTTCG
GGACTCATCACCAATCCCCTTTCAAGAACGTACTAGGCCCCCTCCTTCCTTTTTACGGCGTCATGTCTGC
CTTAGATCAAATTGCCGATCTGTCTCATGCAATCGAGGTCATAGCGAATGAAACCGGTAAGGCTCTCCTG
ATCATCTCCAGTGAACTCGCCTCAGTCAGACTTCTGGCTTTACAGAACAGAGCCGCTTTGGATTTTCTCC
TGGCAGCCCAAGGAGGAACGTGTGCCGTCATCGGTTCTGAATGTTGCACCTTTGTGCCTGAATATAATGC
TACCATTGCTGATATAGTCAATCACCTC

in the original fasta and
>Lophius-piscatorius_GP-Orthoebolavirus-restonense
------------------------------------agagtggctcactccct------------accca
agaagcctcag-------------------ttgaagc-----gggtacggcgagaagattcc--------
--------------gatatcttcgggactcatcaccaatcccctttcaagaacgtactaggccccctcct
tcctttttacggcgtcatgtctgcctta------------------------------gatcaaattgcc
gatctgtctcatgcaatcgaggtcatagcgaatgaaaccggtaaggctctcctgatcatctcc---agtg
aactcgcctcagtcagacttctggctttacagaacagagccgctttggattttctcctggcagcccaagg
aggaacgtgtgccgtcatcggttctgaatgttgcacctttgtgcctgaatataatgctaccattgctgat
atagtcaatcacctc-------------------------------------------------------
----------------------------------------------------------------------
-

in the alignment used here.
I would like to use IQ-TREE -asr to infer the ancestral state of the sequences in the fasta grouped by the branches in the guide tree. If this is possible, how do I get IQ-TREE to match the sequences to the guide tree? If this isn't possible, what other software could I use?
example files

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.

Thanks updated post below.

It's difficult to answer this without example data.
Before we go deeper please try the following.

Do not use '-' in the taxa name, in old phylogenetics algorithms this was banned because it represented an indel. Thus change "-" to "_". I do this all the time automatically.
You are running an ancestral reconstruction with nucleotide sequences.
In this case the first thing to do is run:

     iqtree -s sequence.fa -m MFP -msub nuclear -o example.tre

Find the tree model and the tree and then input that model and the tree into the asr analysis via the -m and then run the exact model via:

    iqtree -s sequence.phy -m GTR+G -asr -te example.tre

-m GTR+G is an example but its whatever output MFP gives you.
Do you see the output tree of the first step is being loaded into the second step to ensure perfect alignment between alignment and tree?
Also why is sequence.phy used when it's fasta format? phy stands for phylip format and those are mix of two different formats.
I personally think this might be about taxa getting confused, but I'd need to see example data with a tree and I'd like to be sure the alignment and the tree perfectly match before a second debugging.

comments
Using phy direct from mafft. I would output the alignment as a fasta and do the whole thing again, unless you dump stuff here for inspection. phylip is the worst for taxa label bugs, the old phylip is terrible - that would absolutely in errors everywhere. The "relaxed phylip format" is still restrictive.
I'd just use fasta, IQ-TREE will take it and it removes one more possibility. Honestly if you dumped your output here this could be cleared up quickly. This approach of try this try that isn't efficient

Post update
Okay the issues are:

The format you are presenting is fasta
The labels on the fasta are radically different to the alignment
e.g.

"> Lophius-piscatorius_GP-Orthoebolavirus-restonense"

> Lophius_piscatorius

The labels must be identical.
There is a complex issue about concatenation for asr but the problem is the alignment in the first instance must be identical to the tree labels.
Once the labels are identical and we have a clear result then complex stuff can begin
There is another problem with the alignment itself, those indels are huge, at a guess this is a 18S style gene, the only gene I know with such massive indels.
The quotation marks MUST be removed. There can be no quotation marks in the alignment file at all.
There is a huge concern that you are not supplying an alignment file to the asr calculation. In that case the calculation will fail.

Its complicated. Have a go with the above, but my suggestion is to ask someone to help you follow the above instructions. This is because you appear to be trying to perform ASR on an unaligned file and that will violate the safety limits of the calculation, in addition the labels are differing between alignment and tree. There is a post-manipulation that might be needed but this is a some way from that stage.
Again unless the above instructions "suddenly click", I recommend you ask for help from someone. The reason is that this is a complex calculation and IQ-TREE is not the simplest way to perform this calculation.

My personal suggestion here is to use Beast (this maybe in Beast2, it's definitely in Beast). I know this is not quite the answer you wanted but there is a point and click option for ASR. You just have to find where the checkbox is within Beauti GUI (I don't remember where it is). The advantage is it will do everything for you in one step, you just need to use a "lognormal" model, under something generic like GTR+Gamma. The alignment format of nexus (or Newick) can be obtained from MAFFT.
With Beast there is no chance of getting a conflict between the alignment file and tree file, Beast does it all at once and the nexus/Newick (file format needed) doesn't have weird issues like "_" versus "-" and there's no alternative formats like Phylip.
It is certainly worth looking and would simplify things.
My concern is - if you get your tree wrong - everything goes wrong no matter how right you've done it.

Thank you. That partially solves it. To answer your question: I'm using an alignment (phy format from mafft), because the sequences are of different lengths. The approach you suggested does work for producing an asr but this one is then sorted by the tree of the sequences rather than grouped by the species tree.

I see @AlicedeSampaioKalkuhl update

Using a fasta instead of a phylip format sadly still results in sequences from one species scattered across the tree in the iqtree -s sequence.fa -m MFP -msub nuclear -o example.tre step. A very clunky solution would be to rerun everything with one sequence per species per run but apart from that I'm not sure.


text:
How to convert SAM file to VCF file. The BCFTools is giving me errors
How do I convert SAM to VCF file.
I tried to convert the SAM file generated using NCBI BLAST to VCF format when I get this error. I do not get any error while using BWA MEM based aligner in SAM format.
Here is the error I get using BCFTOOLS.
bcftools mpileup -Ou --no-reference output.Sorted.bam -o output.Sorted.bcf -d 1
[mpileup] 1 samples in 1 input files
[mpileup] maximum number of reads per input file set to -d 1
Note: The maximum per-sample depth with -d 1 is 1.0x
bcftools: bam2bcf.c:421: bcf_call_glfgen: Assertion `epos>=0 && epos<bca->npos' failed.
Aborted (core dumped)

Please add more context to your question to help us understand what it is that you're trying to do; any additional information you can give about your general study, or what you want to do with the converted file is helpful. In particular, the error suggests there is an issue with the input, so please show the first few lines of your SAM file, including the header information and a few aligned reads.

I tried to convert the SAM file generated using NCBI BLAST to VCF format when I get this error. I do not get any error while using BWA MEM based aligner in SAM format.

Okay, but you haven't provided any more context, you've just explained what you tried to do as an intermediate step (which can be mostly worked out from the command line you already included with the error).

Can you explain what you mean by "convert SAM to VCF"? What exactly are you expecting in the VCF file? Are you trying to do variant calling here? If so, why are you limiting yourself to only one read per position? Why are you using -d 1?


text:
Freebayes giving problem in converting SAM to BCF
$>freebayes -f Subject_1 output.Sorted.bam > output.Sorted.vcf

terminate called after throwing an instance of 'std::out_of_range'
  what():  basic_string::substr: __pos (which is 166) > this->size() (which is 0)
Aborted (core dumped)
```

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

Please add more context to your question. In particular, explain what you want to do, and what you expected to happen. Because your error appears to be due to an input parsing issue, please also include the first few lines of input including header information, e.g. the output of samtools view -h output.Sorted.bam | head -n 50

I tried to use NCBI BLAST output in SAM format as input to convert to BCF format. I converted the SAM to BAM format using SAMTOOLS and then indexed it. The error generated above is while trying to do this. When I try to use FREEBAYES for BWA MEM based aligner in SAM format, it did not generate any error, but the result file was empty. I tried using BCFTOOLS to convert the SAM format file, and it worked without any issues. Am I missing here something for how to use FREEBAYES ?


text:
Are the Gene names in genome annotations for GTFs of C elegans strains given on WormBase consistent?
I am working with doing RNA seq analysis in different strains of C elegans. I wanted to ask if the orthologous genes were given the same name for reference genomes for different strains on WormBase. That would make my life easier as I could compare the expression across strains for the common genes.
By "name" I meant the identifier, like WBGene00067890.
If this is not so, could you please suggest how I can get same gene identifiers for orthologous genes?
Thanks a lot.

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.

It would be helpful to know

what you mean by "name"- like a non-systematic gene symbol ("P53") or a systematic identifier ("NP_000537.1")?
exactly which strains you are working with. Wormbase does not seem to have very many diverse genome assemblies for C. elegans readily available. (Note that I am not a worm/wormbase person, maybe I'm missing something.)

But gene orthology for the main genomes is given directly in TSV files for each genome on wormbase via parasite. My suggestion would be to map those to your expression matrix files and use them to merge datasets using your favorite method in python or R.
Granted, wormbase seems to be going through a gap in funding so some of their services may be down. I had some difficulty clicking around in their comparator genomes.
Alternatively, you might be evaluating expression based on the N2 reference. (Not that this would be ideal.) In this case there would be no issue of orthology, because there would only be one genome/transcriptome assembly.
Failing all of this, you could use one of various tools (here is another option) to annotate any C. elegans assembly with respect to a high quality reference such as N2.

Thanks!! This info is helpful, but let me clarify about the other thing!!

Maybe @eshaan mark the question as "resolved" and ask a new question? That is the etiquette approach on the site.


text:
request for peer reviewed articles using Markov Chains
I am a mathematician setting up an Research Experience for Undergrads (REU) course for next semester. Students have background in modeling with ordinary differential equations (population growth models, predator-prey relations, neuron firing models; stability analysis of linear and nonlinear systems, stable limit cycles). They have expressed interest in learning about Markov Chains, so I'm currently preparing some material on probability and textbook applications of Markov Chains techniques in the life sciences. I'm having more trouble identifying recent-ish peer reviewed important publications where Markov Chain techniques are used in a substantial way for important questions in biology or medicine. Any field of bio/med works, I just want to come up with a few recent (i.e. last decade), impactful examples of how these techniques are relevant.

@M__ if you have a really coherent writeup on BEAST to help OP I do think that would be valuable, given my treatment was extremely shallow and wide rather than narrow and deep. My main concern is that I think it's a mischaracterization to say that MCs are not used much in biology, so I'm trying to make sure that we have the examples to allow OP to explore if they so desire!

@MaximilianPress :-( it's the end of the week and I'm exhausted. I'll see if I regain any energy next week.

There are very many examples of Markov chains in the life sciences. They are foundational models for many applications:

dynamic programming (sequence alignment)
Hidden Markov Models (genomic segmentation, motif finding, or alignment)
population clustering (STRUCTURE)
phylogenetics MCMC (BEAST)
MCMC is rife throughout bioinformatics in general

There are a huge number of publications using Markov Chain based methods in the life sciences, including in recent years, see for example the results of this pubmed search:

It might help to get more information about what you mean when you write about Markov Chains, because clearly 2 different answerers have come up with very different answers! Is there more info on what your students are interested in?


text:
Gene Sequence retrieval (NCBI vs Ensembl)
I wanted to retrieve the sequence for MPL gene (exon 10 on human chromosome 1), which is implicated in chromic myeloid leukemia. The NCBI gene model displays 12 exons for the gene, whereas the ENSEMBL model lists only 4 exons. The literature mentions that the pathognomonic mutations in MPL gene (MPL S505, MPL W515) are found in exon 10. I am able to retrieve the exon 10 sequence from NCBI gene, for which I need to design primers. However, the ENSEMBL DB displays only 4 exons

Which database needs to be used for sequence retrieval if I need to design sequencing primers ?
2.Are the sequences for genes and other parameters different in both DBs (NCBI and ENSEMBL) ?
3.Is MANE (Minimal Annotation of Nucleotide Sequences) useful ?

What gene? In what species? Which two specific genes are you looking at on NCBI vs Ensembl?

Welcome @Bhanani glad you're here from biology SE and hope you get the answer you are seeking

You wrote, "Are the sequences for genes and other parameters different in both DBs (NCBI and ENSEMBL) ?" . NO was my personal understanding @dthorbur differs. MANE is cool BTW. I might put together an answer if I've time. Personally I'd stick with NCBI gene (MANE is a good solution), but the ENSEMBL team will not agree with that and they did have do have a presence here so I'll wait to see if they answer.

@M__ Thank you. I agree...NCBI is user friendly

Please [edit] all this into the question. And make sure to also include the specific NCBI and RefSeq pages you checked so we can know what you looked at. Remember, most human genes have multiple transcripts, and "exon 10" might be included in some, skipped in others, and be the 9th or 11th or Nth exon in yet others.

ENSEMBL gets their data from other sources, including a lot of RefSeq genomes which are hosted on NCBI. I've noticed that when RefSeq publishes updates to some genomes, it can take a long time to be reflected on ENSEMBL, especially for non-model genomes.
You're best bet is to identify which reference genomes are available, and which is most appropriate for your use case. Things like:

Where was the same isolated?
is it annotated?
is it high quality?
what updates are available?
which does ENSEMBL use?

Thank you......

Can you please explain where Ensembl is showing 4 exons for [a transcript of] the MPL gene? Are you looking at the exons for transcript MPL-204 (which is tagged "CDS 5' and 3' incomplete")?
When I look at the transcript annotated in Ensembl as a Ensembl Canonical transcript (i.e. MPL-201), I see 12 exons defined:
https://asia.ensembl.org/Homo_sapiens/Transcript/Exons?db=core;g=ENSG00000117400;mr=1:43331900-43369416;r=1:43336802-43355795;t=ENST00000372470
Sequences / exons can differ between different databases, but I think the primary transcript sequences should be the same (especially if they're both coming from MANE).