text:
Understanding the fasta file format and using edit distance for alignment
I'm a computer scientist teaching algorithms development in the Fall. One of the algorithms we teach is called Edit Distance, and our folklore is that it is used to compare RNA sequences (is this actually true in practice?).
I would like to have students implement the edit distance algorithm and run it on actual SARS-COV-2 sequences, so I'm trying to understand exactly what I get from the GenBank database. I downloaded this one: https://www.ncbi.nlm.nih.gov/nuccore/1798174254
I am looking at the genomic.fna file. So this is apparently the FASTA file format, and the lines that begin with >MN988669.1 ... are comments. I see comments like:
>MN988669.1 Severe acute respiratory syndrome coronavirus 2 isolate 2019-nCoV WHU02, complete genome

Followed by an RNA string. Is this the start of a new sequence for a different coronavirus specimen? So I could have students extract each of these and run edit distance and then produce a dendrogram or something? How do I find more information on where the samples came from? Is this the right file to be using, or should I use the gbff file? And are the PDB files interesting to me at all (I actually do know what PDB files are)?
Also, are there any recommended data sets where we could do something like track mutations in the virus (and see, e.g., that the NYC outbreak originated from Europe and not China)? Are there other useful algorithms / data that might be interesting for students to study in this vein? Specifically interesting to me would be graph-search algorithms, minimum spanning trees, and network flow. Also any NP-complete algorithms that we could run backtracking on. Obviously taking the theoretical study of algorithms to something as currently topical as the coronavirus has pedagogical value.
Thanks
EDIT:
Based on comments below here is what is taking shape.

Have students implement vanilla EditDistance (which there seems to be some disagreement about which algorithms are named what, so let's say insertions and deletions only, which I would call Longest Common Subsequence LCS). Then a variant that also does alignment (i.e. full Levenshtein distance computation, which I would call EditDistance, but Wikipedia calls the Needleman-Wunsch algorithm with a gap penalty of 1), then maybe Needleman-Wunsch with different gap penalties (if someone tells me what would make sense biologically).
Have students implement basic heirarchical clustering / phylogenetic tree generation algorithms a la https://www.ncss.com/wp-content/themes/ncss/pdf/Procedures/NCSS/Hierarchical_Clustering-Dendrograms.pdf.
Have students run their sequence alignment variants and different clustering algorithms on SARS-COV-2 sequences and report on how the choices of parameters in 1 and 2 change the result and therefore potentially the analysis.
Ask some open ended written response questions on what this might mean for society, whether this introduces ethical considerations for an algorithms designer or are they doing just the math, etc.

My learning objectives (as they are now shaping up) are:

Students will understand that just because their algorithm comes with a proof of correctness, it doesn't mean it's the correct algorithm for the job.
Students will understand that different models / parameters to a model result in different outcomes, and thus even computational problems are not purely computational.
Doing theoretical computer science / mathematics is not devoid of ethical considerations.

I would very much appreciate thoughts on the above.

Hi John, if I get a bit of time over the weekend I will answer question. My apologies I can't do this right now.

@Michael Thanks! I look forward to it. I'm realizing how awesome it would be to have a good dataset where students could see, e.g. that the NYC cases originated from Europe, but I don't know enough to get this dataset put together myself, I don't think. Going to poke at it this week though.

I'm assuming what you downloaded is actually a full RNA genome. It is one of many isolates. You could download the rest and do comparisons across them. You can look at also look directly at the RNA (transcription) or convert it into a different space such as amino acid sequence (translation). Here is an example of what that would look like: reddit

Hi John, if you are interested in alignment in general, I would recommend not to take a viral sequence. AFAIK viruses have their own special alignment problems. I think using some eukaryotic gene sequences would do a better job for a class about algorithms.

I think hemoglobin genes are the bio-informatics 101 go-to for such phylogeny/alignment/evolution introductions  :)

i think that in the end is only question of proteins (so what compare should be proteins) –

@KamilSJaron SARS-CoV-2 is fine, it doesn't recombine too much and even then you can align it well.

Your understanding of FASTA format is about right. The type of basic problem you're eluding to we term "sequence alignment"- edit distance might be okay for teaching but in practise we use other algorithms, e.g. you might be interested in the Needleman–Wunsch or Smith–Waterman algorithms. Richard Durbin et al. wrote a great book that covers these any much more https://www.amazon.com/Biological-Sequence-Analysis-Probabilistic-Proteins/dp/0521629713
Tracking mutations etc. requires more than just alignment though, see "phylogenetics" (i.e. building genetic trees) and "variant calling". Also checkout what the nextstrain team is doing https://nextstrain.org/ncov/global
In general, looking for practical applications for your algorithms is great, but be very careful before drawing any real world conclusions about the coronavirus outbreak from such analyses

Thanks! Yes, I certainly wouldn't draw any real conclusions from a basic analysis. I just want students to get a sense of algorithms in the wild, we could even discuss the differences between basic EditDistance and Needleman-Wunsch and Smith-Waterman, and potentially even bring in ethical issues, like asking how choice of algorithm (naive vs. better designed, for example) could have negative societal impacts if used incorrectly...

And I did not know about those, so thanks for the references :)!

I just looked at that book and its first procedure (UPGMA) for building phylogenetic trees looks to be about the same as the hierarchical clustering algorithm that I was thinking of having students implement in this project as well. So here's a quick follow up question: I would assume that the choice of tree construction algorithm can have drastic effects on the output tree. How is it that biologists validate which algorithm to use?

@John Good question re trees. Simulated data might validate an algorithm but it's tricky. A bioinformatian might use a tool that has several models and compare them, e.g. see http://www.iqtree.org/

The correct method is laborious and it is better to give the students a ready made tree at GISAID to investigate the dispersion of COVID-19 into Europe.
However, a quick approach to getting to the alignment and drawing a tree is easy and will readily complement your established teaching method. What this will give you is a very different phylogeny from edit distances and you would explain the matrix differences between the approaches. I suspect NCBI use Jukes Cantor distances.

Go to blast.ncbi.nlm.nih.gov
Select nucleotides
Paste "MN988669" and hit return
This will produce 100 hits
On the menu page select the "draw tree option"
This will produce a reasonable nucleotide tree for 100 COVID-19 sequences
Select "neighbor-joining" rather than "minimum evolution" (this is hierarchical clustering geared to towards mutation rate heterogeneity between taxa)
There are various point and click options including "examine alignment"
You can then select "minimum evolution" and see the changes in the tree (and there are changes) - TEACHING OBJECTIVE 1 & 2
If you want to examine COVID-19 from the European perspective then select a European isolate, e.g. French isolates and input this into the blast. However, European isolates are basically Wuhan original.

I have provided an example below, which I "rerooted" and represents the closest 100 to you sequence using collapsed clade format (your students can undo this, to examine the contents of a "collapsed clade"). The tree shows the dispersion of strains from the Wuhan seafood market.
There's lots of flexibility and the students can do everything inside ½ hour easily and this would complement your approach. The advantage of my approach to teaching phylogeny and regardless of what  you do and how you do it Blast is central to rapidly obtaining alignment data for students and investigators alike. We use different blast options but blasting through is prerequisite to understand the diversity and clean some information on the population structure.

Also, are there any recommended data sets where we could do something like track mutations in the virus (and see, e.g., that the NYC outbreak originated from Europe and not China)?

Yes there are, today this would be Mesquite available here. In my opinion this is a bit advanced. You can track amino acid mutations, which in my view is easier.


text:
Resources to learn genome assembly workflow for small genomes (like viruses)
I have sequencing data of a few samples of a DNA virus.
I'd like to learn de novo assembly of 'short read' data, to construct a scaffold and then count the abundance of each strain in the data.
I heard about SPAdes as a good choice for these kinds of very short genomes. and also BBmap for statistics related to contigs.
I new to de novo assembly could you also recommend introductory material? For example a full workflow post/paper on genome assembly.
I have some experience with RNA-seq and small RNA-seq data.
Note: I am seeking an overview and test data prior reading Spades and BBmap manuals. Thus I need to learn RNA-seq analysis rather than a aligner.
cross-posted in biostars

https://www.biostars.org/p/455717/

@ATpoint why is that a valid reason to close? The internet is not non-redundant

Thanks for your comment. But as you see in the post, I had indicated the cross-posting and its link from the first day that the post was submitted, both here and in biostar page.

The subject area is huge, viral meta genomics, and there are lots of discussion (from a more specific/technical base) on Bioinformatics SE.

I would recommend a specific tool for viral meta genomics rather than using applying a general tool into a specific technical challenge. Here are the ones discussed on Bioinfo SE for de novo metagenomic assembly tested/being applied to viruses. It is not an exclusive list and needs further discussion.
SAVAGE
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5411778/
Baaijens et al. De novo assembly of viral quasispecies using overlap graphs Genome Research 2017 May; 27: 835–848.
S-aligner
https://www.biorxiv.org/content/biorxiv/early/2021/02/02/2021.02.02.429443.full.pdf
Juanjo Bermúdez s-aligner: a greedy algorithm for non-greedy de novo genome assembly 2021 bioRχiv February: 1-7
Rey
http://denovoassembler.sourceforge.net
Sébastien Boisvert et al Genome Biology. 2012 13:R122, December

source: Stack Exchange

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

sort: votes


text:
Explaining the algorithm of RNA folding: what each symbol & value represents?
I need someone to explain to me from the nuts and bolts how algorithms/ maths is used to work out RNA folding. Explain it to me like I am an alien or child.
I am looking at this paper - https://eprint.ncl.ac.uk/240069
All the books on maths, and information theory & mathematical methods in bioinformatics are either too basic or fail to describe adequately what the research paper is doing in a way I can grasp.
I want to understand how ViennaRNA & NUPACK SPECIFICALLY uses maths to describe and work out RNA folding.
What do the math symbols mean, what are they representing, and why? Why do you write it like that?
In the paper I mention above, it describes how it calculates single stranded folding as per below?
a) What does all this mean?
x ∈ {0, 1}
|x |
= 1 if base i is specified in the
design of x to be bound and 0 otherwise. Further, let p
x ∈ [0, 1]|x |
^ I know conceptually it is referring to vectors and components, I know that the idea is to represent RNA folding as a graph structure in order to describe mathematically, but I do not know or understand exactly how or what it is describing in terms of the reality of an RNA folding. Two vertices connected by an edge represent a base-pair etc etc, that I get, but I cant seem to understand any further than that.
Here is the description from the paper:

Single stranded folding. Firstly, we evaluate the ability of a single
DNA strand to fold into its specified target structure (as shown in
Figure 1 top). Using the secondary structure predictor of the
ViennaRNA 2.0 so‰ware suite [13], we calculate the partition function
of all secondary structures that x might fold into.1 For any strand x,
let |x | denote its length and d x ∈ {0, 1} |x | a vector whose
component d x i
= 1 if base i is specified in the design of x to be bound and 0 otherwise. Further, let p x ∈ [0, 1]|x | denote a vector whose i-th
component p x i denotes the Boltzmann probability (obtained from the
partition function) that base i of strand x is paired with another
base. The single-stranded folding score Ssf is then defined as the
normalized Euclidean distance || · || between d x and p x as Ssf(x) =
1 − 1 |x | ||d x − p x ||. (1) Note that 0 ≤ Ssf ≤ 1 and Ssf(x) = 1 if
x folds unambiguously into its target structure.

Do you understand hybridisation and how that is relevant here?

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


text:
Build tree to calculate beta diversity
I am trying to calculate beta diversity* using the UniFrac distance on my data, each sample having OTUs values.
To calculate it, I saw it required building a tree with percentages that represented the OTUs.
However, I was unable to figure what percentage fit each OTU to build the tree from the instructions I found on the internet.
For example, I saw in this link: http://readiab.org/book/0.1.1/3/1#4.1.2,  that they build this tree:

Could you please explain to me how to build the tree?

This is the diversity between local sites and enables the total in a landscape to be calculated (gamma). The tree would therefore appear (presumably) to represent the diversity between sites.

The key thing we need to understand is the nature of the input. If this is Whitikar's beta this likely a species count. It could be something super clever like a 16S metagenomics between sites. I wouldn't have trouble making a tree (its UPGMA) but I'd need to know what the data represented. A, B, C are likely different foci

Can you give a little more information about what tools and data sets you are using? How did you generate the OTU table etc?
If the underlying data is nucleotide data and this is a phylogentic tree. It may well be a population genetic tree based for example in allele frequencies.
For classical phylogenetics then FastTree is a very simple tool to generate a decent phylogeny. If you have the sequences of each OTU, you can align them in MegaX with muscle, or simply using muscle from the command line if you prefer. You can then generate the tree with this command (presuming you are on windows):
FastTree.exe -gtr -nt < aligned_sequences.fasta > otus_tree.tre

I hope that helps!


text:
Finding linkage disequilibrium when I have a kinship matrix but not pedigree
I have a bunch of samples of a non-model organism divided in two groups: parental generation and offspring. I do not have a pedigree so I do not actually know who is a parent of who, and which are the sibships. I have a bunch of SNPs that I can use to create a kinship matrix, but I cannot resolve families just using it.
I want to perform linkage disequilibrium (LD) between the SNPs. I can calculate the correlation between them, but I feel I am losing a lot of power because there is kinship between samples. Is there a program that could take the kinship matrix in account when calculating the LD?
I could, in theory, reconstruct the pedigree from the markers. I run multiple instance of Colony but the results are underwhelming -- most samples are not allocated (admittedly this raises questions on the SNP situation, but this is a whole different issue).

Welcome to the site. Why not just do LD directly?

what do you mean?  I have only a vague idea of actual distances (and surely no genetic distance) because it is a non model organism. All I can do is a correlation

Well given my data -- unphased, no pedigree, it's not I can get much more than a correlation, hence my question -- I can create a kinship matrix, and maybe use it a prior to have some idea of the segregation and the LD, but I have no idea if anyone made some software able to do that

yes, LD is linkage disequilibrium

I originally thought you were working on humans. I spotted "non model organism", that changes things alot, albeit I guess humans are a non-model organism in a way. The question is what the organism. LD has been applied to every kingdom of life, except viruses. Arlequin isn't bad for relatively small to medium data sets (but tricky to get the format right - thats a separate question). If it's a genome sequence then I still suspect its plink. There are a number of recent LD programs BTW.
I did notice Colony was written in Fortran! Wow, the author must be a physicist at ZSL anyway...
If you have a defined number of SNPs - sounds small to medium data set - then I'd use Arlquein http://cmpg.unibe.ch/software/arlequin35/Arl35Methods.html, Arlquin does loads of pop gen BTW its very famous for AMOVA BTW (very cool calculation). HAPLOVIEW has been used for LD, these authors used it here, what surprised me was they avoided using Arlequin for the LD but used HAPLOVIEW. There's a plethora of recent LD algorithms, never used them so can' say, e.g. quickLD, stuff in R.
What I am definitely saying is the LD works on from raw alleles Thats the most important part of this response. I've attempted to describe the rationale below.

If this is human genomes, or large genome data I'd would recommend plink analysis performed directly on your data, as described here https://www.cog-genomics.org/plink/1.9/ld
Using any sort of matrix will remove the linkage disequilibrium signal, IF its a single matrix for all your data. If its multiple matrices you can perform a correlation statistic between the matrices - but its better just to do a straight linkage disequilibrium calculation, e.g. plink.
The reason why correlation matrices will fail is because LD needs to identify incongruence/disparity and whether that incongruence occurs within random expectations.
Basically you are making a phylogenetic tree - thats what the correlation matrix is about, but the tree assumes clonality. LD is about deviation from random expectations (equilibrium), so either a random distribution, e.g. Poison or a randomisation process is needed. What I'm attempting to say is "a tree" (single matrix) and LD are very different calculations. Whilst, technically you can do LD via matrices (plural), its better to do a LD calculation from the raw alleles.

From the comments, it appears a "tree like structure" (ancestors) is wanted to direct an LD calculation.
What you can do is subdivide the populations and perform LD on each subdivision. A tree (here using a correlation or kinship matrix) is a continuum so it's at what point do you call the the subdivision? Its doable but the information is not presented here.
What you can do is use "Structure" - this has been upgraded alot. This is a HWE, I think using Bayes and it designates how many populations (at LD) are present in your data. You can then perform a specific LD directly on each subdivision.
If you compare the output of "Structure",basically Bayesian HWE you should see a correlation with your kinship matrix. What you are inferring is a partial LD parameter that will define ancestry and whilst its doable it has never been calculated in this manner that I'm aware of: its random or it ain't and thats the criteria of significance. This is definitely how LD works in small genomes (everything below higher eukaryotes).

I think there is a misunderstanding -- I have the raw alleles, and I use them for the calculation, the issue being, since I know for a fact I have families in the data (I just do not know the actual pedigree), I could have a much greater power if this information, or a proxy, i.e. a kinship matrix could be included in the information, because to some degree one could theoretically count how the alleles are passed from (putative) parents to (putative) offspring

@user1256923 answered above


text:
Tajima-Nei Distance estimate with BioPerl
I have been trying to estimate Tajima-Nei distance for my data (link below). The output is s phylogenetics distance (pairwise) matrix which accommodates a detailed transition transversion and the Jukes-Cantor correction.
I´m following this protocol from BioPerl: https://metacpan.org/pod/Bio::Align::DNAStatistics
I have 314 sequences in a fasta file and another file with the list of IDs.
Fasta:
>AVP78031.1
----atgttgtttttcttgtttcttcagttcgccttagtaaactc---------------
------------------------------------ccagtgtgttaacttgacaggcag
a----------------accccactcaatcccaattat--actaattcttcacaaagagg
...

IDs:
AVP78031.1
...

I´m using this Perl script to calculate the Tajima-Nei Distance in a pairwise comparison (314 * 314):
use strict;
use warnings;
use Bio::AlignIO;
use Bio::Align::DNAStatistics;

my $file = $ARGV[0];
my $idfile = $ARGV[1];

if ($file eq "" ) {
  $file = "NT_MSA_S_protein.fasta";
} elsif ($idfile eq "" ) {
  $idfile = "NT_ID_S_protein.csv";
}

#### Considerando un archivo
my @contentIDS;

open (LIST, $idfile) or die;
while (my $l = <LIST>) {
  $l =~ s/\n//g; # delete newline
    $l =~ s/\r//g; # delete CR
  next if (length($l) < 1);
  push @contentIDS, $l;
}
close LIST;

#### .... IDs list
my $stats = Bio::Align::DNAStatistics->new();
my $alignin = Bio::AlignIO->new(-format => 'fasta', -file   => $file);  ### $file: MSA file
while (my $aln = $alignin->next_aln) {
  #print "reading...A\n"; ### DIAG
  my $matrix = $stats->distance(-align => $aln, -method => 'Tajima-Nei');
  #print "reading...B\n"; ### DIAG
  ### Obtaining values for each pair (DISTANCE!)
  WL1:
  foreach my $aaa (@contentIDS) { ### ID #1
    WL2:
    foreach my $baa (@contentIDS) { ### ID #2
    next (WL2) if ($aaa eq $baa);
  my $data =  $matrix->get_entry($aaa, $baa);
  #($data = 0) if ($data < 0);
    print "DISTANCE\t$aaa\t$baa\t$data\n";
    } # END WL2
  } # END WL1
}

exit;

#

This script work it when I tried with small data, however, when I try with my real data this is the error message
MSG: Must provide a DNA alignment to Bio::Align::DNAStatistics, you provided a protein
---------------------------------------------------
Can't locate object method "get_entry" via package "0" (perhaps you forgot to load "0"?) at Tajima-Nei_Distance_NV.pl lin$

This is weird because I review my data for ambiguous characters and the characters are in majority "atcg" and on some occasions "n", at least that there are other ambiguous characters (maybe) that represent protein sequence.
I really don´t understand the message because the fasta file is clearly a nucleotide sequence.
Anyone with experience using Bioperl and estimation Tajima-Nei distance?
Link: https://github.com/MauriAndresMU1313/Example_Tajima-Nei_Distance_Bioperl/tree/main
Any comment or help is welcome!!
Thank!

The problem is the code itself.
I suspect there are R (argenine) and Y (tyrosine) in the sequence data. These are common ambiguities but also represent amino acids. The distance estimate has not accommodated these nucleotide ambiguities but is doing a getter setter for amino acids.
It would be easier to perform the same calculation in Mega X than debugging the code.


text:
Differences in webtool and standalone Clustal Omega Percent Identity Matrix
I'm trying to align sequences using MUSCLE in my own computer. To compute Percent Identity Matrix I'm trying to use Clustal2.1 like it's done in https://www.ebi.ac.uk/Tools/msa/muscle/. Even though I'm doing the same thing, I get rounded up values in my computer.
For example, this is the webtool output:
#
#
#  Percent Identity  Matrix - created by Clustal2.1
#
#

     1: JX310755.1  100.00   92.79   65.29   63.98
     2: M97264.1     92.79  100.00   65.66   63.92
     3: AJ292230.1   65.29   65.66  100.00   66.88
     4: KM379144.1   63.98   63.92   66.88  100.00

And this is the output in my computer.
#
#
#  Percent Identity  Matrix - created by Clustal2.1
#
#

     1: JX310755.1     100      93      65      64
     2: M97264.1        93     100      66      64
     3: AJ292230.1      65      66     100      67
     4: KM379144.1      64      64      67     100

I have both installed clustalw using sudo apt install clustalw and by downloading the tarballs from http://www.clustal.org/download/current/clustalw-2.1.tar.gz, with exact results.
The commmand I'm using is clustalw2 alignment -tree -pim.
Any thoughts on why this could be? Thanks in advance.

So, after a year I came out with the solution - as M__ kindly suggested, using Clustal Omega.
First, align usign muscle:
muscle -in fasta -quiet -out alignment
Then, generate percent distance matrix using clustal omega:
clustalo --infile alignment --distmat-out output --percent-id --full

Try clustal omega, instead of clustalw because it is really old. It should report to 2 dp. Variations do occur depending on how indels are treated, but this is just rounding.
You could try
alignment -tree -pim -score

The idea is to convert to absolute rather than percentage and might overcome the rounding error


text:
Parsing PLINK recombination genetic map
I have a Plink genetic map with 4 columns: Chromosome, rs-snp, Map(cM), Position(bp) and approximately 8,000 lines find the recombination rate  [Rate(cM/Mb)] from a HapMapII genetic map of 4 columns: Chromosome, Position(bp), Rate(cM/Mb), Map(cM) ~ 47,000 lines ... some chromosomes will be 6 times larger than this; not all Plink entries will be found in HapMapII, but 99% will.
I need to parse PLINK output and spreadsheet VLookUP is too slow.
Basically I have to find one missing field matching two given fields in a much larger file. I have tried variations of grep on Linux and FindStr in windows but they have not worked.
Thank you!
Here is the Plink map
Chromosome, rs-snp, Map(cM), Position(bp)
19  rs2312724   0.00187 266034
19  rs1020382   0.002681    267039
19  rs12459906  0.010664    276245
19  rs10409452  0.068867    288246
19  rs3813154   0.092406    292082
19  rs895344    0.099705    292901
... 8300 more rows

Here is the HapMapII map
Chromosome  Position(bp)    Rate(cM/Mb) Map(cM)
chr19   253938  0.154179    0.000000
chr19   256859  0.154528    0.000450
chr19   259772  0.154530    0.000900
chr19   260970  0.154697    0.001086
chr19   261033  0.154888    0.001095
chr19   266034  0.347963    0.001870  <-
chr19   266544  1.256622    0.002047
chr19   266888  1.333840    0.002480
chr19   267039  1.542843    0.002681
chr19   267213  2.171166    0.002950
chr19   267614  1.108918    0.003820
chr19   269380  0.711702    0.005779
chr19   276245  1.791471    0.010664
chr19   277717  2.004476    0.013301
chr19   277776  3.133699    0.013420
chr19   279645  4.135431    0.019277
chr19   280299  4.235347    0.021981
chr19   281360  6.164471    0.026475
chr19   287129  6.113884    0.062038
chr19   288246  6.135026    0.068867
chr19   288320  6.136509    0.069321
chr19   292082  8.817796    0.092406
chr19   292583  9.059132    0.096824
chr19   292901  9.062522    0.099705
.... over 46,000 more lines

I need the Rate(cM/Mb) field imported into a version of the Plink map.
Position(bp) 266034 ; and Map(cM) 0.00187 exist in both files; but Rate(cM/Mb) 0.347963 does not. I need to import Rate(cM/Mb) from the larger file into a new version of the smaller file. The columns can be resorted. I don't need the chromosome column. I will add a final column with mutation rate as well.

What you want to do is a join on two columns, these answers show you how to do it. Is that enough for you or I can write a full-length answer later?


text:
PCA analysis and nucleotide diversity pi
I have 25 individuals' whole-genome sequence (collected from 5 different locations) and after doing a variant calling analysis, I found 2734 SNP+INDELS marker (SNPs + INDELS = 2734 ). Then I plotted the PCA plot (pc1 and pc2) using variant calling data, i.e. raw SNPs.
Now I am trying to interpret the population structure.
Here I can see that these individuals from different locations are genetically different and there is no difference between individuals from the same geographic foci.
Are there any other things that I can report?
EDIT: Each circle(except green) represents a cluster of 4 individuals (that overlap with each other) and the green dot represent one individual

Then I calculate nucleotide diversity (π) values (across the whole genome) of each cluster observed in PCA plot: What is best way to show that information?
The mean pi value of each cluster:
1 5.61 E-05
2 4.76E-0.5
3 0.00045
4 2.56e-05
5 4.69e-0.5

Is there any test that should I perform to see if there are any significant differences, such as panmixia, migration or neutrality?
Can someone explain to me how to interpret and report the finding from PCA and π?

What's being shown on the PCA plot? I'd have expected to see the 25 samples, but it looks like you're plotting just the cluster centers, which makes it impossible to evaluate how tightly grouped the samples are. What's the meaning of the inner green circles?

HI, the sample that I plotted overlap each other that's why it is looking like that. Each circle represent (except green) a cluster of 4 individuals (that overlap) and the green dot represent one individual

@NuclearHoagie I made changes to the original post.

when you say you are interested in population structure, what do you mean? why are you estimating pi within clusters? what is the meaning of the clusters? did you place the clusters together or are they the output of some procedure? Is there a reason that you are not using $F_{ST}$ between clusters, if they are pre-defined? Have you considered using STRUCTURE or another tool that does unsupervised clustering (not just PCA projection)?

I'm happy to answer this if @Ashar is still seeking a solution


text:
How can I get properties from each one letter amino acid?
Im using python and wish to get properties from each one letter amino acid?
I mean
M(Met) => Density:1.340, MW:149.21...
The RDkit methods atom.GetSymbol(), atom.GetDegree() permit this and I am looking for an equivalent solution in Python.
Are there any suitable package, for example in Biopython?

Check ProteinAnalysis from Bio.Seq.Utils.ProtParam in Biopython.

RDKit is a Python solution.


text:
Parse RNA variant effect annotations ("r." format)
I've got annotations for splicing variants in a format like this (this is one variant):

Variant: NM_004092.3:c.88+5G>A
Effect: Retention; r.87_88ins1_88+10:p.(Ala31Glufs*23)

I want to extract which acceptor/donor sites are lost or gained (in this case lost causing intron retention). In other words, the genomic location of the sequence that is retained in mRNA.
I parsed the variant HGVS ID using ensembl VEP (like this) without problems, giving me the genomic location of this variant. However, VEP seems to not like "r." annotations for the effect: I tried ensembl VEP again with RNA, protein, and both annotations, none of which worked.
Is there a tool that can parse these annotations for me? I've got too many to do it by hand.

Note on RNA splicing
RNA splicing and acceptor and donor sites:
An introns contains: a 5' donor site, a branch site located near the 3' site and a 3' acceptor site. Both the donor and acceptor  are required for splicing along with the branch site.


text:
Calling for one specific SNP from multiple sequencing runs
Just learning the basics of bioinformatics and bash shell scripting.
Is there a correct way to search for a specific allele with one or even multiple sequencing runs? I am searching for prevalence of the C2416T allele in the SARS-CoV-2 genome within the ORF1a RNA polymerase gene as a mini project.
I downloaded a bunch of SARS-COV-2 sequencing runs, I think I need to convert them to bam file type than I am not sure how to call for one specific SNP, especially in a loop for 5-7 runs. Can this be done with bash?
I was looking into variant calling from a bam file type using vcf and I can't find a way to call for one specific snp on multiple bam files at once.
Ideally I would want to filter the variant out of the .bcf files or use grep since I sort of understand those the most and I’m not working with a large amount of data or a human genome.

Getting SNPs from whole-genome data is usually done via mapping sequencing reads on the reference (bwa-mem or bowtie are two popular mappers) and then using one of variant calling tools (like freebayes, samtools, or GATK). However, I would also look up if there is something virus-specific.
According to this answer, you could get variant calls of a specific region using samtools if you give it a bed file that contains the region of interest (three columns chr    start    stop format). For example using human reference
samtools mpileup -uDl regions.bed -f hg19.fa file1.sorted.bam file2.sorted.bam | bcftools view -bvcg - > RAL_samtools.raw.bcf

You could also specify a bunch of files, not just two and it would do the variant calling on them all.

Answer from @atpoint [transferred from comment]:
Alternatively, you can use samtools view to only retain reads that overlap the region you are interested in and then feed these into the mpileup. IIRC there was something with that -l option in mpileup that caused trouble but I do not remember what it was. I think it still called all variants based on all reads and then simply filtered for the target regions in -l. The view approach would probably be much faster. Just keep that in mind as an alternative.


text:
What is the tP statistic from ANGSD?
We were looking for genetic diversity in the exon level of an MHC transcript and analysing the data under ANGSD (Analysis of next generation Sequencing Data) written in C++
We used ANGSD to calculate Thetas and Tajima D neutrality tests and one of the thetas was a pairwise theta, however I am unsure what this means and how it differs from Tajima's D?
http://www.popgen.dk/angsd/index.php/Thetas,Tajima,Neutrality_tests
We then divided the tP over length of each exon to find genetic diversity in the exon level. However, as I am very interested in statistics, if someone could explain what the tP statistic means exactly it would be highly appreciated.

While I can't explain it to you directly, you might check out the section on theta estimators in Korneliussen et al. (2013) which is cited on the above ANGSD wiki page. The equation for this estimator in this paper traces back to the original Tajima (1989).

I don't know Theta but Tajima D is a population genetics test that assumes constant sample size. When the number of segregating sites (non-conservative sites in an alignment) significantly differs from pi (Nei's p-distance between a pair of sequence) one possibility is selection given a constant population size.
What Tajima found was cool: under neutrality this relationship conforms to a beta distribution. Just to explain that again so the nucleotide differences between every pair of sequences in an alignment (pi) against the total number of possible mutations conforms to a beta distribution for a constant population size under neutrality.
Its a cool idea - like really cool - but the violations - any level of 'extinction' will skew the distribution, population movement and importantly so could a bias in sampling.
So what you get is a pair-wise matrix at 5% critical probability. If your data entirely conformed to Tajima's D it would be an amazing result, but that is not the result you are looking for. So if Tajima D gave a particular pair <0.05  AND it was supported by other tests THEN (sounds like code) it provides a basis to consider selection as an explanation.
Again I don't know Theta but if both tests were congruent for putative selection I'd say that was reasonable. However, Tajima D alone is not a cool result, its cool theory however but not cool in practice. I hope that explains it?
Note: did any one ever follow this up in trees? Cool if they did because there will be a relationship between the balance of a tree (e.g. Bayesian 'Skyline' plots) and Tajima D.


text:
DNA, molecular descriptors
Are there any programs (preferably command line tools) for calculating molecular descriptors of DNA? I am looking for something like Chemopy or RDKit but specifically for DNA.
Thanks in advance!

I would like to perform classification after all, so the more descriptors are the better for me.

Please [edit] your question and add some details. Explain exactly what "descriptors" you are thinking of, tell us what the input to this program would be and also tell us what operating system you are using.

Hi @BalintBiro, thanks for your feedback and answering your own question. Could you kindly mark this as "accepted", it helps with site stats.

Just found the right solution.
iFeatureOmega does exactly what I was looking for. Sequence based and has CLI.
https://ifeatureomega.erc.monash.edu/iFeature2/index.html
Additionally it supports protein, DNA/RNA and also small molecule inputs.

Thanks for your answer on iFeatureOmega, cool. Could you please mark the answer as "accepted"?


text:
How to analyze an IGV alignment
I'm working on a project where I am analyzing the performance of an alignment workflow. My goal is to find regions in the resulting BAM file where there are outstanding discrepancies or anything that indicates my assembly/alignment has "mistakes". I'm working with human genomes from the 1000 genomes project, they are high coverage Illumina paired end reads, each read is 100bp.
My workflow so far:

input paired end FASTQ files (human) into SPAdes to create contigs.fasta file
align paired end FASTQ files to contigs.fasta file with minimap2 to get SAM file
convert SAM to BAM with samtools
sort and index BAM file
load BAM file into IGV and load contigs.fasta as the reference

In the IGV GUI, I have colored alignments by insert size and pair orientation. This has revealed a bunch of insertions, tandem duplications, and inversions.
I am wondering:

Do these colored reads indicate that the assembler made "mistakes" or could they be there just by chance?
Do SNPs indicate anything about how SPAdes/minimap2/samtools performed?
Is there a better way to sort/color the alignments based on what I'm trying to do?

Thank you so much to anyone who has any advice!
EDIT: here is a screenshot of the alignment. For clarification, there are 8 genomes of interest, hence 8 tracks

Background
Integrative genome viewer (IGV) is a well known Broad Institute tool and available here.
The most recent papers are:
James T. Robinson, Helga Thorvaldsdóttir, Aaron M. Wenger, Ahmet Zehir, Jill P. Mesirov; Variant Review with the Integrative Genomics Viewer. Cancer Res 1 November 2017; 77 (21): e31–e34. https://doi.org/10.1158/0008-5472.CAN-17-0337
igv.js: an embeddable JavaScript implementation of the Integrative Genomics Viewer (IGV)
James T. Robinson, Helga Thorvaldsdóttir, Douglass Turner, Jill P. Mesirov
bioRxiv 2020.05.03.075499; doi: https://doi.org/10.1101/2020.05.03.075499

seeing some browser shots would help. In general, a few weird reads is no big deal. A lot of reads showing systematic issues might be more of a problem. It would also help to know what kind of reads. We trust Illumina more for base-level accuracy, and long reads more for contiguity. What species?

It's a human genome from the 1000 genomes project! And thank you, I just added a photo for clarification

thanks- can you share some of the track labels? I assume some of that stuff is gene tracks (e.g. green bars), but there are some coverage plots that don't show the read pileups, which makes me worry that I don't understand what's going on.

perhaps more to the point, there are 3 pertinent pieces of information you can get from this approach: (1) aligning reads against GRCh38, (2) reads against your assemblies, (3) your assemblies against GRCh38 (off the top of my head). I would suggest doing all 3 of these approaches, and comparing them. 1 and 3 can be visualized at the same time, which might be particularly helpful in finding e.g. regions of ambiguity due to repetitive regions. You appear to be doing 2, which might be a good way to identify regions to check by 1 and 3.

What do you mean by track labels? this site makes it seem like it's just the name of the file I'm using. In this case all my files are from de novo assemblies, I only got the FASTQ files off the internet

A filename might be helpful- basically I just don't know what I'm looking at. Even marking it up in MS paint or something would help.

Added! Thank you so much!


text:
Mapping List of Probes/Primers/Short Oligos to a Reference Fasta/q
Could you help me bioinformatics SE people.
So I have been duckduckgo-ing up a storm and even consulted the AI overlords and haven't come up with a solid way to do this, but it is so fundamental to bioinformatics that it has to already exists and I don't need to be re-inventing the wheel...
So imagine you have a fasta/q file with a "reference" sequence:
AATTTTAAAAAAACCCCAAAAAAAAAAAGGGGGAAACCACAAAA
And a list of short oligo sequences:
TTTT, CCCC, GGGG
How would be the best way to map (if that is the right terminology) the short oligos to the reference seq so that a configurable mismatch amount is allowed:
  TTTT       CCCC           GGGG    CCCC
AATTTTAAAAAAACCCCAAAAAAAAAAAGGGGAAAACCACAAAA

Blast, BWA and Bowtie2 from my understanding aren't ideal for this because they aren't built for really short sequences, something about them using seeded reference lookup?
Using Smith-Watterman or Needleman-Wunsch technically gets the local aligning right but the mismach allowance is based on score not X total allowed. Plus, iterating over hundreds of oligos and aligning against dozens of large references seems inefficient.
The closest thing I have found to what I'm trying to do is Sequence Manipulation Suites' Primer Map. This tool is almost perfect except a) it is a web tool, not a module I could import into python or incorporate into a CLI and b) does not appear to offer any tolerance for mismatches.
Am I just searching with the wrong vocabulary? This has to be an already solved bioinformatics process right?

A reference FASTQ is not a thing: A reference sequence is curated in some way, and a curated sequence does not come with quality scores.

Here is the source code of Primer Map You are not restricted to use it as a web tool only. Moreover you can modify it according to your needs.

@RamRS this is why I put reference in quotes lol. Basically the bare sequences are parsed out of a fasta/q file and are used as targets to align all the shorter oligos to.

@M__ I'm trying to determine coverage of a pool of primers against arbitrary target(s)

@haci Thanks! I cloned a copy and am working on dissecting it and figuring out if it has the ability to be modified for mismatching. Still feel like this is something that already exists and I just don't have the right search terms to find it...

@RPINerd I don't quite understand what you mean by "coverage of primers" but you might want to check out primer3. I used this in the past to design/test a large number of primers. You can customize and run via CLI.

@haci taking the example I gave in the OP, the "coverage" would be (4 bp * 4 oligos)/ 44 bp target sequence, so about 36%. And taking start stop positions into account I could also make some statistics on the distance between aligned sites

Both of @haci's solutions are very cool:

Primer Map
Primer3 CLI

I wasn't aware that Primer3 was available as a CLI (command line interface).
The solution I had in mind was in silico PCR one example being Kelpie and the associated paper is here, however I am not exactly sure about the application of the "primer mapping", but most "primer" applications concern PCR.


text:
How to get the number of complete phage genomes available on ncbi?
I am looking to establish the total number of complete phage genomes available on NCBI.
I am not looking for any specific type, but want to understand the total diversity available.
How can I perform the search?

Could you supply more details? What groups of phages? Any in particular? Every phage? The answer depends on what you are looking for.

The idea is to get the total number of complete phage genomes available on ncbi. So yes, I am not looking for any particular kind, just want to know how many are there.

That's easy. Please update your question to clarify this point, because its helps the archives.

Done, thank you for guiding me to make the question better suited for the forum.

Go to the NCBI nuccore database, at https://www.ncbi.nlm.nih.gov/nuccore
Copy and paste the search term:
(bacterial virus[Organism]) AND (3000:110000[sequence length]) NOT bacteria[Organism]

You'll get 74583 results.

Click on the link send to and click the checkbox called file.
It will then download.

Check the downloaded file is in fasta format with ">" to start the ID.

If fasta format is absent, return to the webpage and where it states summary switch that to fasta text

I can't remember if point 4 is needed.
Notes

Sequence length, I am guessing that all phages are between 3000 and 105000 in length, this is to screen out partial sequences and to minimise  bacteria creeping in to the search. I don't quite understand how a phage can be 105 000bp in length, but you can twiddle these parameters to refine the search.
@acvill informs us that giant phages greater than >200 kb have been annoated. Thus a separate search is likely required targeting this genome size.

Anecdotally, (3000:110000[sequence length]) should capture the vast majority of phage genomes, though a growing number of huge phages (>200 kb) are being annotated.

Thank you, that is good expert opinion. I will update the post. I am mystified how a virus can be this size in comparison to their bacterial host genome, but phages are not my thing.

I did not know you could download all the search results at once as fasta.

How about trying this?

Go to https://www.ncbi.nlm.nih.gov/nucleotide. Then type the following command in the search box. It resulted in 23971 items.

(bacterial virus[Organism]) AND (complete genome[ti]) NOT shotgun[ti] NOT plasmid[ti] NOT bacteria[Organism]


text:
use Kallisto in galaxy
I want to use Kallisto for sequence alignment in Galaxy. Its description is:

a program for quantifying abundances of transcripts from bulk and
single-cell RNA-Seq data, or more generally of target sequences using
high-throughput sequencing reads.

However, I found this field empty:

there is no option available for the reference transcriptome.

How I use a reference sequence?
Can I leave the other configuration by default?

this is the link of the tool:  https://usegalaxy.org/root?tool_id=toolshed.g2.bx.psu.edu/repos/iuc/kallisto_quant/kallisto_quant/0.43.1.3

You have to click on the dropdown menu and select "Use transcriptome from history", which will allow you to specify your own file.

I suspect the default will be built-in human transcriptome.
Solution is not tested, so let us know if it actually works.


text:
Steps for trimming aligned peptides using Trimal
I am trimming the aligned peptide sequences and back translating them to codon using Trimal.
Here i am using following cmd's:
trimal -in  OG0008206.fa.aln.out -out OG0008206.fa.aln.out_1st.trim -automated1 -backtrans OG0008206.fa.codon
trimal -out OG0008206.fa.aln.out_2nd.trim -phylip_paml -nogaps

But i am confused about, should i remove the gaps fisrt and then backtransleted the peptides or Backtransleted first and then remove the gaps. Because i am getting different  results (codon) in two attempts. which one i should do fisrt?
Thanks.


text:
Can't load files into PopART
I am trying to create haplotype networks in PopART (Population genetics with Reticulated Trees*) which I have successfully done many times in the past.
Two weeks ago the program was working fine for me but recently I am unable to load files. What happens is that when I go to load a nexus file through the "nex" button, or through the "file/import/alignment" or "file/open" options nothing loads.
There is no error message, instead nothing happens. This occurs with both files I have recently created and with files that I have successfully loaded in the past. I've sent files to other people and everyone has been able to load the files using PopART on their computers. I am running macOS 10.15.6 and have tried turning my computer off and on as well as uninstalling and reinstalling PopART. I also thought that maybe Java wasn't up to date on my computer but I don't think PopArt runs with Java. Any help would be very appreciated!

PopArt builds reticulated HAplotype networks such shown below

Answer from @user438383, converted from comment:
Hi, this seems like quite a niche problem that you are facing, perhaps a bug in the GUI of some kind. Although it's on-topic, I think you're unlikely to find somewhere here who can help (unless the popart devs lurk here). I would suggest emailing the authors, but they point towards the github which mostly seems to be ignored.

Thanks, ill email the authors.

PopArt have a mailing list popartusers+subscribe@googlegroups.com
If you go here you can subscribe to it. This sounds like an issue with Mac compatibility. OS 10.15.6 isn't a version I understand and suspect there's a typo. Currently I believe it is OSX 12.1, thus OS10.12.1
PopArt was written for OSX10.6 and upwards. HOWEVER, there can be issues with major system upgrades so Monterey OSX12, I think was a major revision. PopArt is no longer supported so it would be sensitive to changes in Mac OS.
My suggestion is to try this on a Windows machine, if you get the same error the suspicion is there has been a small change to the input format.


text:
How to convert Sars-CoV-2 ORF1ab codon positions to ORF1b positions?
Some websites mention amino acid mutations as ORF1ab, most however use ORF1a/b instead.
So how can I translate from one to the other? What do I need to subtract from say ORF1ab:4588I to get the corresponding N in ORF1b:<N>I?

ORF1ab:4588I corresponds to ORF1b:187I so one needs to subtract 4401 to go from 1ab to 1b and add 4401 to go the other way.

Just for clarification: 1a or 1b? You reference 1b in the question but answer with 1a.

You're right, I fixed it, thanks!

There is now an Open Source online tool that does the conversion, Codon2Nucleotide by Theo Sanderson.


text:
Python package or CLI tool outputting mutations of sequence with respect to reference
I'm new to bioinformatics, and I'm looking for a tool that takes two FASTA files, one containing a reference, let's call it A.fasta, the other another sequence, let's call it B.fasta, and outputs just the mutations that B contains with respect to A.
I'm open as to the output format. Is there a standard for this?
Background: I'm trying to analyse Sars-CoV-2 sequences and identify mutations with respect to the reference. I've done a google search for a tool that does what I want to but couldn't find something straightaway. Is there no standard tool for this?

I would say parsimony, but thats just me

If it's just a comparison of two sequences that you're looking for, then diffseq from the emboss toolkit should be okay:
https://www.bioinformatics.nl/cgi-bin/emboss/diffseq
I think this tool requires the sequences to be oriented in the same direction. If you discover a lot more variants than you expect, try reverse-complementing one of the sequences.


text:
How to get strain names/ids contained in a multi FASTA file using seqkit?
FASTA files can be very big and unwieldy, especially if lines are at most 80 characters, one can't speed up browsing them by using less with -S to have one sequence every two lines.
How can I extract just the strain names (or sequence names, i.e. the string on the line starting with >) into a list? I'd like to use seqkit since I'm already familiar with it and know some of the flags and idioms it uses.
I googled the question but couldn't find an answer using seqkit - the examples all use awk, grep, sed or similar, which are alright for copy/pasting but not very memorable.

The nice thing about seqkit is that the answer is really simple and memorable:
seqkit seq -n sequences.fasta

This will output a text file containing one sequence name per line:
MPX/PT0001/2022
MPX/PT0002/2022
MPX/PT0003/2022

Thanks for posting back @CorneliusRoemer thats useful.

You can trivially do this with standard *nix tools. Here are a few options:
sed -n 's/^>//p' file.fa
awk 'sub(/^>/,"")' file.fa
grep '^>' file.fa | tr -d '>'
perl -ne 'print if s/^>//' file.fa

The advantage of these is that they will work on any POSIX-compliant system and most other *nix flavors, even if they're not POSIX such as Linux, out of the box without needing to install any external tools.
Another thing to take into account is how fast a command is, especially when dealing with large files. I ran a few tests using the solutions presented here on the human hg38 genome assembly which can be downloaded here and seqkit version 2.2.0 which was downloaded from here. The file was stored locally on my laptop's SSD drive, and the laptop itself is:
$ inxi
CPU: quad core Intel Core i7-6820HQ (-MT MCP-)
speed/min/max: 2896/800/3600 MHz Kernel: 5.17.5-arch1-1 x86_64
Up: 25d 10h 27m Mem: 15035.9/31958.1 MiB (47.0%)
Storage: 476.94 GiB (103.1% used) Procs: 452 Shell: Bash inxi: 3.3.15

To get more accurate results, I used a little function I have which will run the command you give it a specified number of times and then report the average real time across all runs. Here's the code if you want to try it:
timethis(){
    max=$1;
    shift;

    for com
    do
        printf '\nCOMMAND: %s\n' "$com"
    c=0;
    while [[ $c -lt $max ]]
    do
      let c++;
      ( time -p eval "$com" ) 2>&1 | grep -oP 'real.*?\K[\d\.]+'
        done |
        awk -vm=$max '{k+=$1}END{print (k/m)}';
    done
}

And the results. These are the average "real" (i.e. wall clock) times taken across 10 runs:
$ timethis 10 \
  "rg '^>' GCA_000001405.29_GRCh38.p14_genomic.fna | tr -d '>' >/dev/null" \
  "grep '^>' GCA_000001405.29_GRCh38.p14_genomic.fna | tr -d '>' >/dev/null" \
  "sed -n 's/^>//p' GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null" \
  "perl -ne 'print if s/^>//' GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null" \
  "./seqkit seq -n GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null" \
  "awk 'sub(/^>/,\"\")' GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null" \
  "perl -nle 'push @res, $_  while s/^>//; END {print  join \"\n\",@res}' GCA_000001405.29_GRCh38.p14_genomic.fna"

COMMAND: rg '^>' GCA_000001405.29_GRCh38.p14_genomic.fna | tr -d '>' >/dev/null
0.474

COMMAND: grep '^>' GCA_000001405.29_GRCh38.p14_genomic.fna | tr -d '>' >/dev/null
0.583

COMMAND: sed -n 's/^>//p' GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null
5.62

COMMAND: perl -ne 'print if s/^>//' GCA_000001405.29_GRCh38.p14_genomic.fna  >/dev/null
6.408

COMMAND: ./seqkit seq -n GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null
4.326

COMMAND: awk 'sub(/^>/,"")' GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null
13.518

COMMAND: perl -nle 'push @res, 3  while s/^>//; END {print  join "\n",@res}' GCA_000001405.29_GRCh38.p14_genomic.fna
8.455

As you can see above, the grep or rg and tr approaches only took around half a second on average, so that is clearly much faster than any of the other choices. However, seqkit is pretty good since it only took around 4 seconds, faster than any of the other ad hoc solutions.

Seqkit is arguably a bioinfo standard tool ;)

@CorneliusRoemer I guess. I have never used it myself and have been doing bioinformatics for more than 20 years. But by "standard" I mean "standard *nix tools", not bioinformatics tools. I edited to clarify.

@CorneliusRoemer there are weakness on reliance on a developers 'kit'. Flexibility, if you wanted to count the number of pangotype A.1 in SARS-CoV-2, my code is easy to modify to />.A.1./ and introduce $count++ , whereas you would need to search seqkit to see if the developer had foresight. Bugs: are 'kits' completely bug free, how can you be sure? Upgrades: the promised upgrade ain't ever released (very, very common). Distributing your algorithm: this is a headache because you've now got a dependency which you don't have permission to redistribute.

@M__ of course there are use cases for *nix tools. But there are also use cases for *kits. In this case I wanted to use a memorable solution that will not require lookup every time. The two other answers are arguably off topic because they don't address the question - which explicitly asked for seqkit usage. Another *kit would have been more on-topic than *nix tools. I don't mind them, it's just a bit typical for SO, you ask a question and you get a solution that's for a different question ;)

@CorneliusRoemer answers providing different approaches are welcome across the Stack Exchange network (note that this site is not [SO]). Especially when there is already a fine answer with using the requested tool. As for what is memorable, that is entirely personal. You are familiar with seqt so that seems easy and memorable for you. I have been using *nix tools for 20 years, so nothing is more memorable than a simple sed s/foo/bar/. But that's fine! That's the beauty of it: we can all provide answers and people can choose whichever approach works well for them. Everybody wins!

Fair :) though it's seqkit not seqt

Heh, you see how memorable I find it, @CorneliusRoemer! :P

@CorneliusRoemer you've got to agree with terdon: seqkit is far too slow. It is worrying. This would cause problems on the data sets in question where a million strong fasta file is normal. Take your pick but realistically its the end of this discussion, there's no point using an inefficient algorithm against huge data sets, sooner or later they'll be a RAM issue.

That's a good reason to use an alternative tool, I agree, if speed matters. If you use the command line interactively though, terseness can be more important. It depends on the use case I suppose. I'd be curious how fast rg would be. The slowdown in seqkit seems to be related to file access? Large sys time.
Can you add a link on how to download the benchmark data so I can compare on my machine and try rg?

@CorneliusRoemer good catch! The test results I show were using a genome file stored on a network drive, which would of course affect the times. I had assumed that wouldn't be relevant since all tools were accessing the same file over the same network. However, I downloaded a different genome file locally and on my local SSD, seqkit is far more competitive. I'll update in a few minutes with the new test results, but do you have any idea why one tool would be affected by the network I/O times more than another? That's very surprising to me.

Ah interesting! I thought you were the bioinformatician with 20y of experience. I barely have 1 :D
Seqkit could read in chunks? No idea, one should maybe even open this as a bug on Github. By the way there are other similar tools out there.
Please try rg too, it's a Rust port of grep and likely faster than grep. Also best to link to the data file so it's reproducible. Performance could differ on macOS and *nix. Ah yes, that's an advantage of seqkit it could possibly work even on windows...

@CorneliusRoemer Ha! Serves me right for throwing my age around :) I do indeed have 20 years of experience, but I still come from the world of biology, not programming so there will be holes in my knowledge. I have updated the post with the new results, including all data needed to reproduce them, and seqkit now only loses out to grep + tr but is faster than all the others. So much for my years of experience, I have to eat my words now!

@M__ my initial numbers don't hold up when doing the same thing with a fasta file stored on a local SSD. So it looks like seqkit is slower than expected when reading data over a network but is faster than everything except the simple grep + tr when running on a local file.

@CorneliusRoemer if it is chunking it at source, its not a bug. It would make sense. Does it parallelise? Chunking and parallelising would be a technique to speed the operation and minimise RAM issues. Obviously it will fall over via a network.

Thanks @terdon its surely grep/ripgrep. If seqkit is chunking at source, that is really interesting (its clearly not pulling the whole file in and then chunking).

Or a bit more verbose ...
perl -nle 'push @res, $_  while s/^>//; END {print  join "\n",@res}' myfile.fa


text:
Entrez (Biopython) esearch and efetch not returning sequence as expected
I'm trying to use Entrez (through Biopython) to download the sequence of a TMV replicase gene.
I have the following code:
# Lookup ID
search = Entrez.esearch(db='gene', term='Tobacco mosaic virus[Orgn] AND replicase')
read = Entrez.read(search)
idlist = read["IdList"]

# Fetch sequence
search = Entrez.efetch(db='gene', id=idlist[0], retmode='text', rettype='gb')
read = SeqIO.read(search, "genbank")
tmv_replicase = read.seq

But when I run it it throws the exception on the SeqIO.read line: ValueError: No records found in handle.
The IDs provided in idlist are: ['1494081', '1494082'].

Did you try to fetch from the database nucleotide instead of gene ?
# load modules
from Bio import SeqIO
from Bio import Entrez

# Lookup ID
search = Entrez.esearch(db='nucleotide', term='Tobacco mosaic virus[Orgn] AND replicase')
read = Entrez.read(search)
idlist = read["IdList"]

# Fetch sequence
search = Entrez.efetch(db='nucleotide', id=idlist[0], retmode='text', rettype='gb')
read = SeqIO.read(search, "genbank")
tmv_replicase = read.seq

It returns
>>> read
SeqRecord(seq=Seq('GTATTTTTACAACAATTACCAACAACAACAAACAACAAACAACATTACAATTAC...GGG'), id='MN912489.1', name='MN912489', description='Tobacco mosaic virus isolate TMV:BR:SC:02:01, complete genome', dbxrefs=[])
>>> tmv_replicase
Seq('GTATTTTTACAACAATTACCAACAACAACAAACAACAAACAACATTACAATTAC...GGG')

Thank you, this gets me what I wanted. Out of interest though, do you know why using the gene database for this wouldn't work? (I'm new to Entrez).

It is simply 2 differents databases. The ID into nucleotide is not the same into gene. Thus from the website NCBI it is possible to extract a fasta from a gene query it doesn't seem possible with biopython https://biopython.readthedocs.io/en/latest/chapter_entrez.html

If you want to filter your queries from the nucleotide database in order to get only the specific genes you can check my code as an example: https://gitlab.mbb.univ-montp2.fr/bev/target_genes_from_ncbi/-/blob/master/select_query_gene_sequences_ncbi.py


text:
How to get phylogenetic tree from multiple genes?
I constructed a phylogenetic tree using a gene (example - secA). I had to gather the same gene sequence for all the required species from public database-NCBI and then constructed the tree after multiple sequence alignment. I used the MEGA X software.
I now want to build a single phylogenetic tree using concatenated multiple genes, around 20.
How do I go about it.
Do I,

Concatenate all gene sequences for each species and do a multiple alignment or
perform an alignment for each gene and finally concatenate the alignments to get the tree?

Could this be done using MEGA X?


text:
scRNA-seq, 10x cellranger pipelines
I'm starting to do sc-rnaseq using 10x cellranger pipelines, and i add TdTomato sequence to mouse reference genome and add an entry in the gtf. My code is But when I makref, it remindered that:
Writing genes GTF file into reference folder...
Invalid number of columns in GTF line
Can anyone tell me the reason?

Do you have a fasta and gtf file for your custom gene? It sounds like you can just add them below your genome.fa and genes.gtf files respectively with cat command.

I have add  fasta and gtf file to  genome.fa and genes.gtf files, but cellranger did"t work ,and I update my problem and related pictures,can you heip me find the reason?

The error messages you shared within your question and within what you have pasted is different. The first error is about line 1794438, which is 20 . exon 1 1431 . - . gene_id "TdTomato"; transcript_id "TdTomato"; gene_name "TdTomato". As you can see, the first field of that line is not the same as your gene name that was used in the fasta file, but is "20". I guess it is also a good idea two double check the spacing between different fields in your GTF, which is causing the problem and should be tab-separeted.

yes,thank you.  I should check the spacing between different fields in  GTF.but the first field of 1794438 line is 20,( 20 . exon 1 1431 . - . gene_id "TdTomato"; transcript_id "TdTomato"; )why you say is "20", I did't find.

thank you all for your reply and help, the reason for the bug is when I copy gtf code from windows to linux, the tab changes,and gtf was not tab-separeted.Thank you again!

@sophia,I see you replaced your initial question and made it into a new one, Please leave the original question intact, and start a new question for your other problem.

Please don't edit the question so that it asks something different! That invalidates the existing answers. If you need help with the step, post a new question. And please never post screenshots of text like that! We need to be able to copy your data and code so showing us an image is essentially useless.

As was suggested by @benn, you will need to add the TdTomato sequence to the end of your fasta formatted genome file. If your TdTomato sequence is in fasta format already, cat will do the trick. If not, just go the end of the file, start a new line and add a header line starting with ">" then immediately the name of your gene, no whitespcaes in between, something like ">TdTomato". Then add your sequence in a new line.
The format of the GTF file is explained in the image you posted, the trick is to use the exact same name (in column 1) that is used in the genome file you will create at the previous step (no ">" though). Again, this line, in which the different fields (columns) will be separated by tabs, will be added to the end of your GTF file.

Thank you for your detailed answer,  and I have add  fasta and gtf file to  genome.fa and genes.gtf files, but cellranger did't work ,and I update my problem and related pictures,can you heip me find the reason?

thank you all for your reply and help, the reason for the bug is when I copy gtf code from windows to linux, the tab changes,and gtf was not tab-separeted.Thank you again!


text:
Which correlation method to compute the correlation score between different clusters of Sc-RNAseq data?
I'm analyzing single cell rna-seq data and trying to compute the correlation score between different clusters. Wondering how to choose the correlation method("pearson" (default), "kendall", or "spearman") of the Cor function?
I'm a new comer in bioinformatics and know little about these arithmetics. Can anyone give me some advice or show me the code to decide which method to choose for my data matrix?
Below is my gene expression matrix: rows are gene names,and columns are clusters.I use the default Peasrson method to do analysis,but I really don't know if it is right.

CorOb.cor.exp <- as.data.frame(cor(CorOb.av.exp))
CorOb.cor.exp$x <- rownames(CorOb.cor.exp)
neupeptide.cor.df <- tidyr::gather(data = CorOb.cor.exp, y, correlation,CorOb.cor.exp$x)
a <- ggplot(neupeptide.cor.df, aes(x, y, fill = correlation)) + geom_raster()+scale_fill_gradientn(colors = c("blue", "white", "red"))+
  theme(axis.text.x = element_text(angle = 45, vjust = 1,
                                   size = 12, hjust = 1))+
  coord_fixed()

a

Cross-posted: https://github.com/satijalab/seurat/issues/2198

I don't think github issues are the best place to ask for general help such as this.

This question does not have a definite answer. Rather, the answer depends heavily on the assumptions you wish to make regarding your data.
Given a pair of random variables $x$ and $y$, Pearson's rho has a value of
1 when there is a linear relationship between $x$; i.e. $y = x\beta + c$, where $c$ is a constant, and there is no noise.
Note that the slope, $\beta$, and intercept, $c$, of this relationship
do not affect the correlation coefficient. With Spearman's rho, we first compute the rank
of each element in $x$ and $y$. We then calculate the Pearson correlation between
these ranks.
Kendall's tau correlation coefficient also deals with the ranks of $x$ and $y$,
but is slightly more complex. The most relevant difference in this case is likely
that tau is better in the case of ties, and may be more robust for low sample sizes.
In most cases, Kendall's and Spearman's are likely to give very similar results
If you wish to assess only linear relationships between expression profiles,
Pearson's is the most appropriate. If you wish to allow for non-linear effects,
Spearman or Kendall is most appropriate (Kendall if the data contains ties).

Disclaimer: This is not a complete answer but more a comment, posting as an answer because of character restrictions for comments. Also, I am not a statistican but a molecular biologist who did quite some reading/courses on stats.
Two scRNA-seq tools, namely SingleR and RCA (Reference Component Analysis, Li et. al, 2017) compute correlations of bulk transcriptomic data to single cells; for this they use Spearman correlation and Pearson correlation respectively (The two tools use the same bulk transcriptome reference). I am using these two tools for different purposes so I did not really compare their results, however, single cell clusters that I compute using the reference components (calcuated with RCA) is in good agreement with single cells labels computed with RCA (adjusted rand index of 0.81).
In addition to that, please take a look at scMatch: a single-cell gene expression profile annotation tool using reference datasets, they claim (scMatch computes bulk - single cell correlations like SingleR and RCA):

From this comparison we found Spearman’s correlation coefficient
  outperformed Pearson’s correlation for all six libraries.

Pearson correlation evaluates two continuous variables. Spearman correlation is often used to evaluate ordinal variables. Kendall correlation is a nonparametric test for two continuous variables.


text:
bowtie2 options when mapping stranded single end reads
I'm trying to map RNA-Seq reads generated using the NEB Ultra Directional kit is first strand reversed protocal. I'm using bowtie-2.3.5_1 to align against the mus musculuous (mm10) genome.  The reads are single ended.
It seems that --fr and --rf are relevant only for paired end options. I tried these options as well as not specifying either (default).
E.g.
$ bowtie2 --fr \
        --threads 20 \
        --very-sensitive \
        -x /reference/mus_musculus/mm10/ucsc_assmebly/illumina_download/Sequence/Bowtie2Index/genome \
        -U trimmed_out/my_trimmed.fq.gz \
        -S bowtie_fr_1.sam

$ samtools sort -n -o bowtie_fr_sort_1.sam bowtie_def_1.sam

I then compared the output of bowtie_rf_sort.sam bowtie_default_sort.sam and bowtie_fr_sort.sam.  They are basically identical excluding the header.
Question:
Does strandedness matter in this alignment? Perhaps I should care about the --nofw/--norc options? Am I missing something here?

Are you taking transcripts as the reference? --norc might be what you are looking for.

I'm using mm10 as my reference (downloaded from Illumina).

Then you need an RNA-seq mapper, not bowtie2. STAR and HISAT2 should have options to enable strand specific mapping. I am not familiar with those, though.

Could you please elaborate how bowtie is not an rna-seq mapper? Per their website "Bowtie 2 is often the first step in pipelines for comparative genomics, including for variation calling, ChIP-seq, RNA-seq,"

Bowtie2 is not aware of splicing. For typical RNA-seq applications, you will want to use a splice-aware mapper, such as star and hisat2, which is specifically designed for RNA-seq. Tophat2 does use bowtie2 for mapping, but it is invoking bowtie2 in a non-standard way and is generally thought to be superseded by star and hisat2 anyway.

I repeat what user172818 says...don't waste your time on bowtie or tophat.  Use STAR or HISAT2, or a pseudoaligner like kallisto with a transcriptome.


text:
Which candida albicans fasta and gff file should I use for alignment?
I downloaded the gff and fasta files for candida albicans from http://www.candidagenome.org.  I want to use hisat2 to align some fastq files against it.
There are many files in each download :
gff :
A22_ForcheSNPs.gff
A22_Historic_Assemblies.gff
A22_Jones_PMID_15123810_Polymorphisms.gff
A22_Jones_PMID_15123810_Polymorphisms.vcf
A22_Unannotated_transcripts_Bruno_et_al.gff
A22_Unannotated_transcripts_Sellam_et_al.gff
A22_Unannotated_transcripts_Tuch_et_al.gff
C_albicans_SC5314_A22_current_features.gff
C_albicans_SC5314_A22_current_features.gtf
C_albicans_SC5314_A22_current_features_with_chromosome_sequences.gff.gz
C_albicans_SC5314_A22_current_intergenic.gff
C_albicans_SC5314_haplotype_variations.gff
C_albicans_SC5314_version_A22-s07-m01-r179_features.gff
C_albicans_SC5314_version_A22-s07-m01-r179_features.gtf
C_albicans_SC5314_version_A22-s07-m01-r179_features_with_chromosome_sequences.gff.gz
C_albicans_SC5314_version_A22-s07-m01-r179_intergenic.gff

fasta :
C_albicans_SC5314_A22_current_chromosomes.fasta.gz
C_albicans_SC5314_A22_current_default_coding.fasta.gz
C_albicans_SC5314_A22_current_default_genomic.fasta.gz
C_albicans_SC5314_A22_current_default_protein.fasta.gz
C_albicans_SC5314_A22_current_not_feature.fasta.gz
C_albicans_SC5314_A22_current_orf_coding.fasta.gz
C_albicans_SC5314_A22_current_orf_genomic_1000.fasta.gz
C_albicans_SC5314_A22_current_orf_genomic.fasta.gz
C_albicans_SC5314_A22_current_orf_plus_intergenic.fasta.gz
C_albicans_SC5314_A22_current_orf_trans_all.fasta.gz
C_albicans_SC5314_A22_current_other_features_genomic_1000.fasta.gz
C_albicans_SC5314_A22_current_other_features_genomic.fasta.gz
C_albicans_SC5314_A22_current_other_features_no_introns.fasta.gz
C_albicans_SC5314_A22_current_other_features_plus_intergenic.fasta.gz
C_albicans_SC5314_version_A22-s07-m01-r179_chromosomes.fasta.gz
C_albicans_SC5314_version_A22-s07-m01-r179_default_coding.fasta.gz
C_albicans_SC5314_version_A22-s07-m01-r179_default_genomic.fasta
C_albicans_SC5314_version_A22-s07-m01-r179_default_protein.fasta.gz
C_albicans_SC5314_version_A22-s07-m01-r179_not_feature.fasta.gz
C_albicans_SC5314_version_A22-s07-m01-r179_orf_coding.fasta.gz
C_albicans_SC5314_version_A22-s07-m01-r179_orf_genomic_1000.fasta.gz
C_albicans_SC5314_version_A22-s07-m01-r179_orf_genomic.fasta.gz
C_albicans_SC5314_version_A22-s07-m01-r179_orf_plus_intergenic.fasta
C_albicans_SC5314_version_A22-s07-m01-r179_orf_trans_all.fasta.gz
C_albicans_SC5314_version_A22-s07-m01-r179_other_features_genomic_1000.fasta.gz
C_albicans_SC5314_version_A22-s07-m01-r179_other_features_genomic.fasta.gz
C_albicans_SC5314_version_A22-s07-m01-r179_other_features_no_introns.fasta.gz
C_albicans_SC5314_version_A22-s07-m01-r179_other_features_plus_intergenic.fasta.gz

I'm pretty sure I'd like to use the version *A22-s07-m01-r179* since it seems like it is a version that will is not transient like the *current* files.
QUESTION :
I am looking at bulk gene expression derived from an RNA-Seq experiment run on an Illumina sequencer.  Which gff / fasta file should I use for alignment with hisat?

What kind of sequencing did you do? Targeted? A gene panel? An exome? WGS?

Just gene expression using rna-seq on an Illumina sequencer.  Grinding up cells and getting the bulk gene expression

OK, please add that to your question, but make sure to actually explain what you are looking at. The file you need will be different if you're looking at all genes or specific genes. You can always use the main file with all chromosomes, but if you know you only care about 5 genes and all or on chr2, there is no point in using the whole thing.

First of all, you almost certainly want the current files. The convention in pretty much all sequence repositories is that foo.current is always the link pointing to the most recent and most up to date data. Those are what you want and not some "stable" version, you want the most accurate and that will usually be the current. In any case, my bet is that you will find that the *A22-s07-m01-r179* files are exactly the same as the corresponding *current* files since the current ones will just be symlinks.
Next, in order to align your data, you need the genome. The genome is this one:
C_albicans_SC5314_A22_current_chromosomes.fasta.gz
That is the file that has the chromosomes in fasta format, so that's the one you want for alignment. Now, if you have done RNAseq and only care about annotated, known genes, you could speed things up a little by using one of the more targeted files (just scroll down on the page you linked to for an explanation of what each file has), but using the whole genome won't be wrong. And at 8.4M compressed, that's a tiny genome so I don't think it's worth trying to get much more speed out of what will already be a pretty quick process.
As for which GFF, that will again depend on what you want, but most likely you want either of these, depending on whether you want GFF or GTF:
C_albicans_SC5314_A22_current_features.gff
C_albicans_SC5314_A22_current_features.gtf

The main message is that you usually want the "current" files, just make a note of what version the "current" files are pointing to if you need to report that.

The refseq is Candida albicans SC5314.
I assume you are performing a fasta reference based assembly.
Its 8 chromosomes are  NC_032089.1 to NC_032096.1 inclusively from chromosome 1 to chromosome 7 (NC_032095.1) to chromosome R.
Its here.
Most of the files you downloaded are SC5314. So I dunno it depends what the project is re-annotation, reassembly?
The vcf is a very interesting file, particularly if you are sequencing your own variants of Candida.
Whats in all the files:

C_albicans_SC5314_A22_current_orf_coding.fasta.gz proteins ignoring introns - useful
C_albicans_SC5314_A22_current_chromosomes.fasta.gz is going to be NCBI files I described.
C_albicans_SC5314_A22_current_default_protein.fasta.gz orfs splicing out the introns

... list goes on really it depends on the project. I would recommend current only. A22-s07-m01-r179 will likely be the previous assembly / annotation.
I guess the single file you want is C_albicans_SC5314_A22_current_chromosomes.fasta.gz for referenced based assemblies (presumably of Candida isolates).


text:
Using multidimensional scaling to visualize protein sequences by functionality
I have a multiple sequence alignment (MSA) of protein sequences, with which I have performed multidimensional scaling to visualize their clustering.  Using R, I used the seqinr package to create an identity matrix (using dist.alignment).  Is there a similar method that is more "generalizable"?  Specifically, I'm interested in classifying the amino acids by some functional characteristics, such as polarity, instead of having the alignment list the amino acids by their identity (A, G, I, etc.).
Here's a very short set of fake input data, with columns 2-4 representing amino acid residue positions:
a <- c("Strain 1", "Strain 2", "Strain 3", "Strain 4")
b <- c("Polar_hydrox", "Nonpolar", "Nonpolar", "Polar_hydrox")
c <- c("Nonpolar", "Nonpolar", "Nonpolar", "Nonpolar")
d <- c("Nonpolar_cyclic", "Nonpolar_cyclic", "Polar_amide", "Polar_amide")
e <- c("Polar_amide", "Acidic", "Polar_amide", "Acidic")
x <- data.frame(a, b, c, d, e)

Is what I'm asking a feasible/correct use of MDS?  Or should I be looking into a different method entirely?  I attempted to use Gower's formula to create a dissimilarity matrix on the original dataset (i.e. using the actual amino acid identities instead of the functional groupings), but compared to what I got using dist.alignment, the MDS plots look different.
Edit: Gower's formula and the dist.alignment function are obviously calculating the distance matrices differently, but I don't know the details on the different methods.  Would using Gower's formula be a good approach for calculating distance matrices for the example dataset?
Edit: clarified functions used for calculating distance matrix.

If you use different distances metrics it is expected to get different MDS plots. What do you mean by more generalizable? What is your biological goal?

Thanks for your response.  My goal is to see whether sequences appear more related to one another when considering what the functionality of their amino composition is, versus the raw amino acid composition.  I'm interested in the phenomena that not every kind of amino acid change will affect immune escape equally.

The question is then why do you settle for Gower's formula and not other distance metrics like the Manhattan distance? It is as good as any other distance metric as far as I can tell.

Some of the first MSA analysis, eg of codon usage of bacteria were performed using matrix factorization (see old papers from Des Higgins). His group also published discriminative methods for finding functional residues using supervised correspondence analysis

Wallace IM1, Higgins DG BMC Bioinformatics. 2007 Apr 23;8:135.
  Supervised multivariate analysis of sequence groups to identify
  specificity determining residues.

The method is available in the R package ade4 and Bioconductor package made4
Also have a look at the phylogenetic R package adephylo and several related packages

I personally don't think this will work well for some variables, primarily because you are not mapping conformational epitopes, where amino acids will cluster in physical space by their charge. In particular, the charge is critical to retaining this conformational structure.
However, hydrophobic/hydrophilic characters are far more amenable to linear style epitopes that this analysis is primed towards.
I don't recommend a multi-dimensional array for this analysis, because you have loads of information, but rather a 'dictionary' (Python), hash (Perl) (or 'object' [any]) or XML, which retains all the information for a given sequence and then strips this in preparation for your multivariate statistic.


text:
Prevent Snakemake from changing wildcards: sample for every rule
EDIT: I posted this as a bug, perhaps wait for the result of that? github
All,
I am using Snakemake to build a pipeline. To define input and output names, I constructed them based on a filename base (from the samples table), a read indicator (like R1, R2) and the fastq file extension, like so:
samples = pd.read_csv('samples.tsv', sep='\t').set_index('samples', drop=False)
reads = ['R1', 'R2']
fastq_file_extension = _001.fastq.gz

rule fastqc:
    input:
        os.path.join(config['Directories']['fastq_raw_dir'], '{sample}{read}' + fastq_file_extension)
    output:
        html=os.path.join(config['Directories']['fastqc_dir'], '{sample}{read}_fastqc.html'),
        zip=os.path.join(config['Directories']['fastqc_dir'], '{sample}{read}_fastqc.zip')
    wrapper:
        "0.49.0/bio/fastqc"

However, during pipeline execution already I can see strange things happening, for this rule it reports:
[Fri Feb 21 12:39:57 2020]
rule fastqc:
input: fastq_raw/0053_P2017BB3S20R_S2_R1_001.fastq.gz
output: qc/fastqc/0053_P2017BB3S20R_S2_R1_fastqc.html, qc/fastqc/0053_P2017BB3S20R_S2_R1_fastqc.zip
jobid: 10
wildcards: sample=0053_P2017BB3S20R_S2_R, read=1

As you can see, Snakemake pulled the R part of the read to the sample and set the read to 1 (should be R1). But this still works because it produces output files as Snakemake expects.
But the next rule, trim_galore_pe:
rule trim_galore_pe:
    input:
        [os.path.join(config['Directories']['fastq_raw_dir'], '{sample}' + read + fastq_file_extension) for read in reads]
    output:
        fastq_read1 = os.path.join(config['Directories']['fastq_trimmed_dir'], '{sample}.1_val_1.fq.gz'),
        fastq_read2 = os.path.join(config['Directories']['fastq_trimmed_dir'], '{sample}.2_val_2.fq.gz'),
        report_read1 = os.path.join(config['Directories']['trim_galore_qc_dir'], '{sample}.1.fastq.gz_trimming_report.txt'),
        report_read2 = os.path.join(config['Directories']['trim_galore_qc_dir'], '{sample}.2.fastq.gz_trimming_report.txt')
    params:
        extra="--retain_unpaired"
    wrapper:
        "0.49.0/bio/trim_galore/pe"

This goes wrong, during pipeline execution Snakemake reports:
[Fri Feb 21 12:40:37 2020]
Error in rule trim_galore_pe:
jobid: 7
output: fastq_trimmed/0053_P2017BB3S20R_S2_.1_val_1.fq.gz, fastq_trimmed/0053_P2017BB3S20R_S2_.2_val_2.fq.gz, qc/trim_galore/0053_P2017BB3S20R_S2_.1.fastq.gz_trimming_report.txt, qc/trim_galore/0053_P2017BB3S20R_S2_.2.fastq.gz_trimming_report..txt
conda-env: /home/nlv24077/projects/snaqs/snakemake_envs/1b94430b
cluster_jobid: 11904662.osios
Error executing rule trim_galore_pe on cluster (jobid: 7, external: 11904662.osios, jobscript: /home/nlv24077/temp/RNASeqPairedEnd/.snakemake/tmp.ooypuklc/snakejob.trim_galore_pe.7.sh).
For error details see the cluster log and the log files of the involved rule(s).

This is understandeble because the file fastq_trimmed/0053_P2017BB3S20R_S2_.1_val_1.fq.gz does not exist, it should be fastq_trimmed/0053_P2017BB3S20R_S2.1_val_1.fq.gz, somehow Snakemake added an extra underscore, which you can see in the execution output, under "wildcards" (0053_P2017BB3S20R_S2_).
[Fri Feb 21 12:39:57 2020]
rule trim_galore_pe:
input: fastq_raw/0053_P2017BB3S20R_S2_R1_001.fastq.gz, fastq_raw/0053_P2017BB3S20R_S2_R2_001.fastq.gz
output: fastq_trimmed/0053_P2017BB3S20R_S2_.1_val_1.fq.gz, fastq_trimmed/0053_P2017BB3S20R_S2_.2_val_2.fq.gz, qc/trim_galore/0053_P2017BB3S20R_S2_.1.fastq.gz_trimming_report.txt, qc/trim_galore/0053_P2017BB3S20R_S2_.2.fastq.gz_trimming_report..txt
jobid: 7
wildcards: sample=0053_P2017BB3S20R_S2_

So how can I prevent Snakemake from changing the samples and read wildcard in seemingly arbitrary ways?

I think you want to prevent snakemake from interpreting your wildcard strings as regular expressions. For this, escape each wildcard value using wildcard_constraints. E.g., towards the top of your Snakefile add:
wildcard_constraints:
    read= '|'.join([re.escape(x) for x in reads]),
    sample= '|'.join([re.escape(x) for x in samples])

(Personally, I use the wildcard_constraints directive as above quite liberally to constrain every wildcard)

Does this mess with how you can use rules.fastqc.output.html as the input for other rules? It gives me a lot of errors for a downstream step (multiqc).

As far as I know, this doesn't cause problems. If you don't show the error it's difficult to tell what's going on.

Placing two wildcards directly behind each other will make it impossible for snakemake to determine where one wildcard ends and the other begins. So it is simply taking the longest possible match for the first wildcard (sample) and then one character for the second (read). This will also cause the _ to be included in the wildcard value sample.
I'm assuming that this same pattern is used for the inputs of another rule to refer to the trim_galore_pe output. Hence why the _ is included in the wildcards there.
You might be able to fix this, simply by adding a _ in-between {sample} and {readgroup} in your first rule (and any other rules using this same pattern):
rule fastqc:
    input:
        os.path.join(config['Directories']['fastq_raw_dir'], '{sample}_{read}' + fastq_file_extension)
    output:
        html=os.path.join(config['Directories']['fastqc_dir'], '{sample}{read}_fastqc.html'),
        zip=os.path.join(config['Directories']['fastqc_dir'], '{sample}{read}_fastqc.zip')
    wrapper:
        "0.49.0/bio/fastqc"

I did this but Snakemake conjures up a new _ out of nowhere:
Building DAG of jobs...
MissingInputException in line 2 of /rst1/temp//RNASeqPairedEnd/rules/trim_galore_pe.smk:
Missing input files for rule trim_galore_pe:
fastq_raw/0053_P2017BB3S19R_S1__R1_001.fastq.gz
fastq_raw/0053_P2017BB3S19R_S1__R2_001.fastq.gz

But why is is impossible for Snakemake to separate them? They are not in {} for nothing right? I don't understand why Snakemake would want to change my wildcards from anything it gets from the samples table. What would be the use? It seems to me that this is just a bug? Or am I really fundamentally misunderstanding Snakemake?

Regarding my first comment, one might think that Snakemake now added and extra _ so I can remove the existing one perhaps? But no, it then removes both:
Building DAG of jobs...
MissingInputException in line 2 of /rst1/temp/RNASeqPairedEnd/rules/trim_galore_pe.smk:
Missing input files for rule trim_galore_pe:
fastq_raw/0053_P2017BB3S19R_S1R2_001.fastq.gz
fastq_raw/0053_P2017BB3S19R_S1R1_001.fastq.gz
So it's either 0 or 2, not 1, as one would expect.

So I guess I am misunderstanding this completely, When I change the input of trim_galore_pe, it complains of missing files, but when I remove the input entirely it complains that the fastqc rule is missing input files, even though I did not touch that rule. Is there any way to more clearly see where errors originate exactly? I feel completely lost now.

@Freek Could you share the entire Snakefile? It's quite hard to see what's going on exactly without it.

Ok, I think I have found the reason, I was using code as well to create the samples_table from the same config and it was adding and removing the _ there as well. So this method does work, my apologies. It remains unhandy that I can't directly use wildcards together, But I could perhaps define a separator in the config file to solve this and keep flexibility..


text:
Combing multiple SRA Runs into a single Sample, BioSample or Library
On GEO/SRA some datasets have multiple Runs per Sample, for example this one (96 SRR numbers for 12 BioSample/GEO_Accession/Library Name/Sample Name), and the authors provide their RNA-Seq counts for 12 samples/condition see here.
I know that an experiment involving Next Generation Sequencing sometimes requires that more reads are generated than fit on an NGS chip lane, or even more than fit on a single Next Generation Sequencer. So a sequencing library can be split over multiple lanes or multiple machines. This leads to multiple sets of raw data which all belong to the same library. For example to reach the right sensitivity in transcript abundance measurements, you need all the reads of the entire library analyzed together. At least, that is how I understand it. A Run is a set of reads, a Library may consist of multiple Runs, right? Here is some more relevant discussion on Reddit, and here is another one. Looks like you can merge but must be careful, there can be (quality) differences between the Runs. Moreover, NCBI themselves say here that they simply sum the counts of multiple runs pertaining to a single GEO Sample (RNA-Seq). It also looks like in the example dataset the number of bases is only sufficient for normal RNA-Seq analyses of all Runs of a Sample are combined together (you want 25 million reads or so for polyA based RNA-Seq).
I have two questions about this:

I think that I can just download the SRA files (prefetch), turn them into fastq files (fasterq-dump) and then sort them and cat them together, right? Or is there a better/recommended way?

What constitutes a single experiment? Is it the BioSample, GEO_Accession, Library Name or the Sample Name? I'm thinking GEO_Accession, since that is what is used on GEO to designate a "sample", right? I need a reliable indicator for what happened to single sample...

Looking forward to hearing how you all deal with this, highest regards.

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.

1, yes, it's ok to merge fastq, but don't sort the combined fastqs, because sometime later on if you need to run pseudo-aligners (eg. salmon), it requires random read id.
A more common way is to treat them separately, (in qc, alignment, filtering, and quantification). Then make sure they look same or close to identical in PCA, then you are safe to merge the final counts to increase the read depth (essentially increase power).
2, GEO samples, refer to the accession with GSMxxx or SAMNxxx, it can contain multiple runs (SRR). If you use SRA Run Selector, and put your PRJNA of interest, you can see the structure of the dataset.

Also, I found that Salmon has options to use multiple libraries as input, see here: https://salmon.readthedocs.io/en/latest/salmon.html#providing-multiple-read-files-to-salmon
I think that is the option I will use.

The rule of thumb is: do not merge FASTQ files. Strange things like same sample different read lengths in runs or different quality encodings while thankfully rare will mess up any downstream analysis.
Also: no point of duplicating input data, creating giant single files and in turn slowing down the analysis.
Process the FASTQ files separately, make sure that runs are from the same sample/same read length, then merge sorted BAMs.

Salmon does not produce BAM files, but I guess I could simply add up all the counts from the runs, then re-normalize. Allthough Salmon can handle multiple fastq-files per library so there is that option...