# Question
Title: Export multiple variables from shell script into a Nextflow output channel
Within a Nextflow workflow, some process runs a shell script that creates multiple variables, e.g.

```
var_a=$(echo "1")
var_b=$(echo "2")
var_c=$(echo "3")
var_d=$(echo "4")

```

These variables contain statistics (e.g. number of input reads, quality scores, ...) that I want to keep until the end of the pipeline, to then save them and perform QC. So I'd like to export them out of this process and carry them to the end of the workflow. To be clear, the number of variables is not unknown, it's just too many for convenience if I need to list them all every time.

I understand that I can use the `env` qualifier or an intermediate file to export the variables individually:

```
process get_vars {
  output:
    tuple env(var_a), env(var_b), env(var_c), env(var_d)

  shell:
  '''
  var_a=$(echo "1")
  var_b=$(echo "2")
  var_c=$(echo "3")
  var_d=$(echo "4")
  '''
}

workflow{
  get_vars | view
}

```

with output:

```
[ea/a145b4] process > get_vars [100%] 1 of 1 ✔
[1, 2, 3, 4]

```

However that becomes very impractical if I have a lot of such variables (that I do want to continue passing along to later processes, along with a much longer tuple of other relevant variables).

My question: is there a convenient way to "encapsulate" all these variables in the output channel, for example with a Groovy map similar to the `meta` map in nf-core? So I could directly write something like:

```
process get_vars {
  output:
    val(my_vars)

  script:
  '''
  var_a=$(echo "1")
  var_b=$(echo "2")
  var_c=$(echo "3")
  var_d=$(echo "4")
  '''

  my_vars = [:]
  my_vars.a = var_a
  my_vars.b = var_b
  my_vars.c = var_c
  my_vars.d = var_d
}

workflow{
  get_vars | view
}

```

which obviously fails because `var_a` is not defined outside of the `'''` block.

# Answer
I am not sure if i get what you mean by encapsulate, you can give your tuple a name using emit?

```
process get_vars {
  output:
    tuple env(var_a), env(var_b), env(var_c), env(var_d), emit: vars

  script:
  '''
  var_a=$(echo "1")
  var_b=$(echo "2")
  var_c=$(echo "3")
  var_d=$(echo "4")
  '''
}

workflow {
  get_vars()
  get_vars.out.vars.view()
}

```

That gives:

```
 [1, 2, 3, 4]

```
> 2 votes

---
Tags: nextflow
---


# Question
Title: Why does DNABert use overlapping k-mers as input?
I was reading through the DNABERT paper and found that their input features were k-mers. This is essentially using rolling window features.

Why train a transformer-like model on sequences of the form ${a\_n,..,a\_{n+k}}$ rather than training it on the complete sequence ${a\_1,...,a\_{m}}$ with positional encoding?

To directly quote the paper:

> Instead of regarding each base as a single token, we tokenized a DNA sequence with the k-mer representation, an approach that has been widely used in analyzing DNA sequences. The k-mer representation incorporates richer contextual information for each deoxynucleotide base by concatenating it with its following ones. The concatenation of them is called a k-mer.

There doesn't seem to be any ablation study in this paper showing that this really is of any benefit.

I would appreciate any more information: I just don't see how the 3-mers {ATG, TGG, GGC, GCT} contain more information than the sequence {A,T,G,G,C,T}, 1-mers, with positional encoding.

# Answer
> 1 votes
NLP and vectorisation is used to trap non-linearity which is highly applicable to sentences written in English. Vectorisation of k-mers for an overlapping sliding window is weird because it's trying to assess the functionality of non-coding sequence, i.e. "silent" and non-cds mutations will be the dominant mutation class, but it would identify localised co-mutations, i.e. within a k-mer.

The only application I can think is stuff like identifying promoters or mutations in promoters, i.e. non-coding regions with a known function and can be trained against appropriate independent experimental observations, e.g. gene expression. If the k-mers were large it might work for RNA secondary structure because it would trap compensatory mutations. Vectorisation for individual nucleotides, i.e. dumbie variables, would not be successful at trapping compensatory mutations.

In summary it has niche application.

BERT is a famous large language model BTW.

---

If k-mer =3 that's about non-linearity of proteins, but it depends how the sliding window is deployed. It depends on the step size of the sliding window.

To try and answer this better, proteins fold back on themself triggering "co-evolution" or more accurately "covariance" between non-linear sites. This defines the triplet codon:

* in machine learning this is needed because it will struggle to find it by itself
* it will help it identify co-evolving mutations at non-linear locations across the gene

It might also identify selection events better.

**Summary** If you've sat down and thought the equation the RNN is attempting to unravel its mind-boggling. More worrisome does it really do it? It's really complicated to simply discuss, but when you give it key information the ML/DL gets there quicker. The authors will be aware of that in their training and testing.

Honestly, we're at the limits of reasonable application of ML/DL without going into RasNet type infrastructure. My system is different the idea is the same - give the training model the key biological information, it gets to the answer quicker, a lot quicker.

# Answer
> 0 votes
\[I'm not an expert in this field, this is only my "best try". Maybe you could get better mathematically motivated answers on ai.stackexchange\]

First note that, as pointed in the DNA BERT paper, k-mers are a very common approach to working with sequences. This is true both for older approaches as well as newer Transformer-type models. For example the recent Nucleotide Transformer, DNAGPT or MuLan-Methyl also use k-mers.

While I'm not aware of a published systematic comparison between k-mers and single bases, I'm sure I've seen papers that vary `k` as a hyperparameter (so would notice if `k=1` gives best results), and I have to believe someone somewhere at some point has thought about not using k-mers.

There are some Deep Learning papers that use individual one hot-encoded nucleotides, for example SpliceAI, but they typically only look at relatively short, fixed-length, sequences (*relatively*: this paper does go up to 10k).

As to the "why", my (limited) understanding is that there are two factors. If you think of NLP, a token is typically a word from the language. Taking each nucleotide as a token means you consider a language with only 4 words, that is a very limited amount of information. It does make sense to use more meaningful tokens: to respond to your comment, even ViT uses patches (and does not attempt to work on individual pixels). Based on our knowledge of biology, we can expect that most sequences have meaningful "motifs" that are made of a few consecutive nucleotides. If you feed the nucleotides to the Transformer, it has to use attention to detect these motifs as well as their long-range interaction. If instead you use a "patch" (or k-mer) that captures entire motifs, this can be learned by an initial linear projection and you only use the attention for long-range interactions.

Finally, the other reason is computing power, as hinted to in the DNAGPT paper:

> In the non-overlapped k-mers strategy, the shift is equal to K, resulting in N/k tokens for an N-length sequence and improving the efficiency by k times.

Maybe taking 1-mers could work just as well as 6-mers, but the number of relationships to consider increases even more, a non-overlapping k-mer strategy does decrease the of interactions to consider, and an overlapping k-mer strategy where each patch is preprocessed (e.g. with a linear model) can also increase efficiency (so to speak the model does not need to consider short-range interactions which are already accounted for).

---
Tags: k-mer, deep-learning
---


# Question
Title: How to find amino acid sequence of a given protein
Is there a way to look up the amino acid sequence of a given protein? For example, what amino acids are used to produce Amylase? I spent hours googling but to no avail. If that's impossible, how can I find the sequence of A, C, G, U of the mRNA that instructs the ribosome to produce Amylase?

# Answer
Usually you go and search in a database. NCBI and Uniprot (specialized on proteins) are ways to go here.

These two offer tons of other related information, links to other databases, advanced functionalities like similarity search based on sequence (e.g. blastp). etc.

To narrow your search, you would need to specify which species you are interested. Here a quick result for amylase in humans, using NCBI and Uniprot. The nucleotide sequence is linked in both cases.
> 2 votes

---
Tags: proteins, sequence-analysis
---


# Question
Title: Genetic relationships between Cyanobacteria: terrestrial vs aquatic
I believe one of the site members is an expert in Cyanobacteria hence the question.

I'm looking for a comprehensive 'phylogenetic tree' of the phylum Cyanobacteria and want to identify all the **terrestrial** members/groups. For example, *Nodosilinea* is one of the terrestrial members.

Part of my rationale for the question is that the Orders can be polyphyletic. From an outsider's point of view this is a difficult taxonomic group to understand.

Orders

> Chroococcales Chroococcidiopsidales Gloeobacterales Nostocales Oscillatoriales Pleurocapsales Spirulinales Synechococcales Incertae sedis

**Background** Cyanobacteria is a massive phylum and important to both terrestrial and aquatic habitats, but what I need to understand is how terrestrial members are distributed across the tree. The reason they are interesting (both aquatic and terrestrial) is that they are efficient at removing CO2 from the atmosphere. Reducing a carbon footprint is an area of growing commercial interest.

# Answer
Ah! Finally, a question I can answer.

There has been some pretty good work about this recently. In particular, I would reccomend this paper by Chen et al. , which groups cyanobacteria into terrestrial / freshwater / marine and...well, the title says the rest. Some are edge cases by necessity (e.g. some coastal biofilm-forming cyanobacteria may be more similar to terrestrial ones than to many marine ones) but that's just nature.

Chen, MY., Teng, WK., Zhao, L. et al. Comparative genomics reveals insights into cyanobacterial evolution and habitat adaptation. ISME J 15, 211–227 (2021). https://doi.org/10.1038/s41396-020-00775-z (https://rdcu.be/dmJpm)

I regularly use this in my work, mainly because they have their trees posted on iTOL (in the 'Data availability' section) and this makes it really easy to search for specific strains.

Which brings us the the main problem. Right now there are about three main ways of separating the cyanobacterial phylum into orders. There is the NCBI database classification (largely based on 16S rRNA), like for example the Chen et al. paper uses. This may sometimes appear polyphyletic when overimposed on a genome tree. Then there is the GTDB classification (https://gtdb.ecogenomic.org/), which is based on genomes. This can be seen in Alcorta et al. 2020 (Alcorta J, Alarcón-Schumacher T, Salgado O and Díez B (2020) Taxonomic Novelty and Distinctive Genomic Features of Hot Spring Cyanobacteria. Front. Genet. 11:568223. doi: 10.3389/fgene.2020.568223) - on average, the tree shape is similar to the one in Chen et al., but the orders make more sense. The GTDB classification is good phylogenetically, and well-integrated with the database (for example, you can just get the list of all genomes/metagenomes available for an order). But because it's not focused on the taxonomy (the naming of things), the actual names may not always fit with the Bacterial/Botanical Code. One research group has been kind enough to go through the absolute mess that is cyanobacteria taxonomy, and now this paper exists.

Strunecký, O., Ivanova, A. P., & Mareš, J. (2023). An updated classification of cyanobacterial orders and families based on phylogenomic and polyphasic analysis. Journal of Phycology, 59(1), 12-51. (https://onlinelibrary.wiley.com/share/A5UXIHBUCIH9JRTIWFSE?target=10.1111/jpy.13304)

As far as I've seen, their orders more or less match with GTDB, just under different names (except for GTDB Cyanobacteriales with is split into more orders in this paper). They also have a much larger 16S rrNA phylogeny constrained by the genome one.

Hower, that solves only part of the problem. Because at the smaller scale (genus, species), you still have strains labeled by their morphologies, and culture collections are (for understandable reasons) hesitant to change their naming system. For particularly visually non-descript groups like 'Synechococcus' or 'Leptolyngbya', there is a good chance that a particular strain could be nearly anywhere in the phylogenetic tree. I would always get a 16S rRNA (if not already avaiable), BLAST it, and then check where related strains appear in a genome tree.

Edit: I have also heard of Cydrasil (https://www.cydrasil.org/), an online tool which can place a 16S rrNA sequence within the context of the cyanobacterial tree, but I haven't tested it yet.

Hope that helps!
> 2 votes

---
Tags: phylogenetics, phylogeny, bacteria, taxonomy
---


# Question
Title: Remote NCBI's Blast Perl API: maximising hits returned via command line
**Question and Background** I'm using NCBI's Blast Perl API to interrogate Genbank for example `blastp`. I don't get the throughput using their point and click web stuff.

I tried Biopython's `Bio.Blast import NCBIWWW` but its under access restrictions from the host, notably attempting to use it during the day (works sporadically). The API however gets through and runs as fast as a web search ... that's good.

Here's an example of the command I'm running:

```
CMD=Put&PROGRAM=blastp&DATABASE=nr&HITLIST_SIZE=2500&ALIGNMENTS=2500&QUERY=%3Emytest%0AMSILANKDTR%0A

```

The above query gives default Blast output. What I'm trying to do is increase the number of hit recovered to the maximum (5000). The options for this are `HITLIST_SIZE`, `ALIGNMENTS`, `DESCRIPTIONS`. In the above example I'm attempting a **hit\_size** of 2500 ... but the same output is returned as for default. More generally I can't get any of the "number of hits/descriptions/alignments returned" options to work, the default is returned each time.

Has anybody solved this?

Typical options for the API are given below

**Typical options**

```
Parameter   Definition  Type    CMD Allowed values
QUERY       Search query    String  Put *   Accession, GI, or FASTA.
DATABASE    BLAST database  String  Put *   Database from appendix 2 or one uploaded to blastdb_custom (see appendix 4)
PROGRAM     BLAST program   String  Put *   One of blastn, blastp, blastx, tblastn, tblastx. To enable megablast, use PROGRAM=blastn&MEGABLAST=on.
FILTER      Low complexity filtering    String  Put F to disable. T or L to enable. Prepend “m” for mask at lookup (e.g., mL)
FORMAT_TYPE Report type String  Put, Get    HTML, Text, XML, XML2, JSON2, or Tabular. HTML is the default.
EXPECT      Expect value    Double  Put Number greater than zero.
HITLIST_SIZE    Number of databases sequences to keep   Integer Put,Get Integer greater than zero.
DESCRIPTIONS    Number of descriptions to print (applies to HTML and Text)  Integer Put,Get Integer greater than zero.
ALIGNMENTS      Number of alignments to print (applies to HTML and Text)    Integer Put,Get Integer greater than zero.

```

---

Just to feedback on the answer and comment

* Using `-remote` by locally installing blast as recommended by @Scott. This is a valid answer so definitely upvoted and recommended its further upvoted. I had tried this but not mentioned it in the question (really sorry), but it times out due to traffic issues. The API overcomes traffic issues. Again, @Scott's is a valid solution (considering I'd not explained I'd tried it).
* More details ...

The Perl API is here https://blast.ncbi.nlm.nih.gov/docs/web\_blast.pl

Line 95 below ...

```
$args = "CMD=Put&PROGRAM=$program&DATABASE=$database&QUERY=" .
$encoded_query;

```

This is changed to,

```
$args = "CMD=Put&PROGRAM=$program&DATABASE=$database&HITLIST_SIZE=$maxhits&QUERY=" . $encoded_query;

```

Also need to add an additional line at line 73

```
$maxhits = shift;

```

From the command line,

```
perl ./web-blast.pl blastp nr 2500 ./example.fa

```

(Note, The variable order is switched in the program).

example.fa - e.g.

```
>rleguminosarum
WGFDGSSTQQAEGRSSDCVLKPVAIYPDPARTNGALVMCEVMMPDGVTPHASNARATILDDEDAWFGFEQ
EYFFYQNGRPLGFPEQGYPAPQGPYYTGVGYSNVGDVAREIVEEHLDLCLAAGINHEGINAEVAKGQWEF
QIFGKGSKKAADQIWMARYLLQRLTEKYGIDI

```

I'm proficient in Perl5, i.e. there's no error in changes to the code in the API, albeit my expertise is object oriented and high performance Python.

---

**Limit on the request** ... I'm testing that now.

First results *eeek* @terdon might be right ... :-(

# Answer
While I haven't used your exact method above, using the BLAST command line tools return the desired number of target sequences

```
blastp  -query query.fa \
        -task blastp-short \
        -db nr \
        -max_target_seqs 5000 \
        -outfmt 6 \
        -out query.out \
        -remote

```

query.fa

```
>mytest
MSILANKDTR

```

The resulting `query.out` file does contain 5000 unique hits.
> 2 votes

# Answer
Right here's the response from NCBI,

> Please send your future blast related questions to blast-help@ncbi.nlm.nih.gov directly. The gb-admin alias is for submitters with questions related to their sequence deposition/submission.

> For your issue, you should look at the following two parameters: DESCRIPTIONS Number of descriptions to print (applies to HTML and Text) Integer Put,Get Integer greater than zero. ALIGNMENTS Number of alignments to print (applies to HTML and Text) Integer Put,Get Integer greater than zero.

> Do please note that asking for this level of hits saved/returned could be problematic for the search. They could make the search exceeds the cpu time allocation and lead to error.

I tried this ... it was the same as "HITLIST\_SIZE".

I figured the answer ...

**Why its happening?** The API uses a default parameter for the **maximum number of hits*. It's imposing a limit on the number of hits returned, but this wasn't stated in the documentation of the parameters.

**Whats the limit?** It's \>70 but \<= 100. I suspect the limit is 100 hits. Thus anything above this its set at 100.

**Summary** The reason why the API is working and other approaches simply time out is NCBI is imposing an upper limit on the number of hits for the API and the Blast servers are not permitting requests above this threshold at peak time. Thus because other approaches are not imposing the restrictions on the NCBI server, the API gets through, whilst the other methods (blastp directly and Biopython) are timing out.

**NCBI are recommending the API is because it is automating a 100 hit cut-off and the servers are deprioritising limits beyond this**

I have filed a note of advice to NCBI because it should be clear in the documentation there's an upper limit on the HITLIST\_SIZE parameter. It took me at least a full half day to figure this out (which isn't cool for limited documentation).

---

**Proof of concept (untested)** Thus, using @Scott's solution

```
blastp  -query query.fa \
        -task blastp-short \
        -db nr \
        -max_target_seqs 100 \
        -outfmt 6 \
        -out query.out \
        -remote

```

This will work with the `-max_target_seqs 100` and thats all the API is doing.
> 0 votes

---
Tags: blast, perl, api, parameters
---


# Question
Title: I am trying to create a subset of 10k variants from 25-30 unmapped contigs of a g.vcf file including the header
My objective is to take a g.vcf.gz file and from 25-30 unmapped contigs with titles like "NW\_020192317.1", I want to make a subset of ~10k variants from each of the unmapped contigs and make one final g.vcf file that includes the header from the original g.vcf.gz file.

From the post here, I have tried running:

```
samtools view -bo output.g.vcf -s 123.4 inputfile.g.vcf.gz NW_020192317.1

```

and I keep receiving the output:

```
Aborted (core dumped)

```

I'm not too picky about the number of variants per contig, as long as it is in the thousands, that's why I tried using ".4"

Alternatively, I tried using the code:

```
samtools view input.g.vcf.gz "NW_020192317.1:1-2" > output.g.vcf

```

and I recieved the same response of "Aborted (core dumped)" with nothing else. I am confused what is wrong here and why it is throwing an error.

# Answer
> 2 votes
After some guidance from my advisor, this command allowed me to get the output file I needed.

```
bcftools view -o newvcfFileName.g.vcf -s 123.2 locationOfTheFile chrom1 chrom2 ... chrom23(mito DNA) --force-samples

```

I listed every chromosome's title until I typed --force-samples. I found it simpler to copy and paste the names in.

---
Tags: vcf, samtools, gatk, subset
---


# Question
Title: Is a classification tree appropriate method to use for my project?
I have about ~3000 genomes of organisms and presence (recorded as 1)/absence (recorded as 0) data of ~300 genes in them. Among 3000 genomes, we have experimental data on 20 organisms, which were represented by 2 different phenotype groups (a, b). According to a previous study, the 300 genes I mentioned are believed to be associated with phenotypes 'a' and 'b'.

To predict the phenotypes of 2980 genomes with no phenotype information, I did hierarchical clustering of 3000 genomes based on the presence/absence of the 300 genes. I cut the resulting dendrogram at a certain height that allowed distinction of clusters that included either group 'a' or 'b', but not both. This gave me 120 clusters, including clusters with no group 'a' or 'b' organisms, clusters with at least one group 'a' organisms, and clusters with at least one group 'b' organisms. Among these 120 clusters, there are 15 clusters that had either group 'a' or 'b' organisms.

I am questioning which collections of genes are responsible for causing phenotypes 'a' or 'b'. To address this, I am thinking about assigning different names for the 15 clusters that contain group 'a' or 'b' organisms, and make a classification tree using these 15 clusters. Does this sound like an appropriate application of classification tree approach? If not, could there be a more sound way to address this question?

Thanks!

# Answer
> 1 votes
This isn't a cool strategy because dendrograms are very poor representations of the relationships between genomes, for example they use UPGMA style metrics fail to account for mutation rate variation between genomes and prone to long-branch attraction artefacts.

I would recommend doing whole genome alignments using e.g. muscle5 -super option and then performing a formal phylogenetic tree using e.g. iqtree or even just MEGAX (distance methods) and using bootstrapping to identify the robustness of clades comprising a and b.

If two clades are supported by bootstraps \> 80% then your hypothesis looks fine. At present your downstream analysis could be an artefact of the dendrogram phylogeny.

---
Tags: clustering, machine-learning
---


title: efetch multiple records at once

text:
# Question
Title: efetch multiple records at once
If I want to efetch two records, I can do the following:

```
efetch -db nuccore -id AB610939,AB610940 -format fasta

```

I have hundreds of records to fetch. Can I put all the IDs in a file and efetch. If so, how?

# Answer
Yes, `efetch` can be given a file with IDs. See `efetch -h`:

```
$ efetch -h | grep -m1 file
  -input         Read identifier(s) from file instead of stdin

```

So you can pass your file of IDs (e.g. `ids`) directly:

```
$ cat ids
AB610939
AB610940

$ efetch -db nuccore -input ids -format fasta
>AB610939.1 Synthetic construct RNA, spike-in microarray control RNA CRM 6204-a, clone: 6204-a-500-1
CTCGACTAGTTAATACGGTACAGGATAACCGATCGGCTTGCAACATAACGGCGTTAAGAATGCGGGAGTG
CAGTTTCCGATTCTCACATCAATCGCCAATAAGGCCTTGTCGCAATATAGACTCAACGGTTCTAGTAGCT
GATCGGTATTACGTGACGCAACCGATTAGACATGCACAATTCCTTGGTCGCTATACTACGGAAATCGTCA
GGTACTATAACCCGTCGCAGGCCTAATACGTGTCGTCACATCGCCAACCTATCGTCAGTCGGAAAGACGT
TGCTGTCTACCATCGAAACTATTTACCGCTCCGAGATTCACGAGTACGAACTCACGAGGAAGTTGCCCTA
TGTAAGGTATCACTCCAGGTACTGCGCCGATAGTACCAGGTGATCAAACGGTTGCAAGAAGGCCACGACG
TATCGGGCTCTTTAGACGTACGCTCGAGATTAAACGCGCACTGATTCACTTTAGCCCGGAATGTCTCGGT
GCGATGTAGA
>AB610940.1 Synthetic construct RNA, spike-in microarray control RNA CRM 6204-a, clone: CRM6204-a-500-2
AGACTAAATCTCGGCGTCGGTTCATACGCGCGATCGTTTGCTGTCAGGGCATACTCGAATCCGGACTCCG
ACAATTATAGGCCATCCTGAATAGCCGATCATGCGAGTCACGATAAGGCAGGCTCTGCGATATCCCGATA
TACTGGAGAAGCTGAATCCCACCTAGAGCGAACTGTCAGAGGATCGACCTCAGGCTCGCTATCATCATAA
CGGCGGACGACCTGTGTCACATTCCGAACGCTACGTGACGATATTATCTGTCGAAAGGCATAGAACGCCG
GTCAATATCCTGCGGCATTCTCTTTATCACCGGCTATAACTACTAGGTTCCGCAGATATAGACTGCGCAC
GGAACATGTAGATAGATCGAGTAGGGTAGCGATTTAACGACTCGACTTACAGACAGAGACGTAGAACGTC
AGACGAGTGGTATGCCCACCAGAGGCGATACAGGCTGTACCTGCGTAGCACTAGAGTCGTGCGTCATGCG
GACCCTATCT

```

If that were not possible, another option would be to put all your IDs as a comma separated string. If your IDs are in a file named `ids`, one per line, you can convert that to a comma-separated string like this:

```
$ tr '\n' , < ids | sed 's/,$/\n/'
AB610939,AB610940

```

That assumes GNU `sed`. If your `sed` doesn't handle `\n`, you can use this instead:

```
printf "$(tr '\n' , < ids | sed 's/,$//')" "$(cat ids)"

```

Armed with this, you can use command substitution to pass the comma-separated list to `efetch`:

```
efetch -db nuccore -id $(tr '\n' , < ids | sed 's/,$/\n/') -format fasta

```

If there are limits to te number of IDs you can fetch at once, or as may well be the case, if this results in a command that is too long, then you can use a loop as suggested in another answer.

---

However, be aware that while looping seems to work for `efetch`, for reasons that are not clear to me, the `esearch` developers have chosen to make it consume stdin. This means that the obvious, straightforward approach fails with `esearch`. See https://unix.stackexchange.com/q/682748/22222 for details.
> 3 votes

# Answer
If you're happy looping through each id in turn and each id in the single file (`idfile.txt`) is separate by a newline:

```
while read -r x; do
  efetch -db nuccore -id  "$x" -format fasta
done < idfile.txt > fetchout.fa

```

*Speculation* If the ids are not separated by a newline but a comma on a single line that might give the simultaneous efetch you're wanting. I strongly suspect that each id would need to be surround in quotation marks in such a file. Again this is just speculation and the above solution is preferred IMO.
> 2 votes

---
Tags: edirect
---


# Question
Title: Using average of sequence similarities when delimiting genera - how to deal with outliers?
Sometimes I am trying to see where a (bacterial) genus ends and where another one begins based on 16S rRNA phylogenies. Of course, the length and support of the branch matter a lot, but this question in particular is about sequence similarity.

As far as I understand it, in general, sequences within the same genus are expected to have roughly \>94-95% sequence similarity at the 16S rRNA level. But that leaves out a lot of variation. For example, in doodle A below, there may be \<90% similarity between sequence 'a' and sequence 'e' ; however, there is always \>95% similarity between, say, a and b, b and c, c and d, d and e...So all would belong to the same genus.  Previously I had considered just comparing every sequence with every other sequence in my group of interest, get an average and a standard deviation, and see if the average \> 94-95%. However, it occured to me that if I have a tree shaped like that in doodle B (a couple of closely-related sequences from one genus and one very distantly-related sequence from another genus), I could get the very same average (and only slight changes in the standard deviation).

It feels like there should be a simple mathematical property to distinguish between A and B based on similarity values, and/or point out outliers like sequence 'a' in B, but I do not know what that would be. Help would be appreciated!

Mock dataset A:

| sequence1 | sequence2 | sequence similarity (%) |
| --- | --- | --- |
| a | b | 99 |
| a | c | 93 |
| a | d | 92 |
| a | e | 95 |
| b | c | 99 |
| b | d | 98 |
| b | e | 97 |
| c | d | 95 |
| c | e | 94 |
| d | e | 88 |

Mock dataset B:

| sequence1 | sequence2 | sequence similarity (%) |
| --- | --- | --- |
| a | b | 90 |
| a | c | 90 |
| a | d | 90 |
| a | e | 90 |
| b | c | 99 |
| b | d | 99 |
| b | e | 99 |
| c | d | 99 |
| c | e | 97 |
| d | e | 97 |

# Answer
> 1 votes
The basic question is where do genera start and stop within phylogenetics and the 16S 5% rule? This is not a trivial question within molecular taxonomy, especially for bacteria.

1. The key is to define the tree structure for 16S
2. Define the tree structure for other genes

**Background and tree stats** The difference between Tree A and Tree B is referred to as 'balance'. Tree A is unbalanced, whilst Tree B is balanced and that is a formal description. I would personally use the term Tree A was "heavily paraphyletic", or "extensive paraphyly" - at a technical level that is accepted, I have used frequently. Equally I'd describe Tree B as a "singular monophyly". The statistical basis for tree balance is called the Colless value.

**"Colless statistic"**. It's fairly easy to calculate. The reason it is not often used is because the tree needs to be phylogenetically robust (below). Not much importance is given to it because if a branch isn't robust thats a headache for this statistic (it can't tell the difference between robust and non-robust).

**Difference between branches** The statistical difference between taxa A and taxa B is measured by bootstrapping (I'm sure you know this), where \>= 80% is bifurcation is robust.

Colless statistic might be used, but its fairly obvious, bootstrapping is needed and **iqtree** is best package for this operation.

---

The second more difficult part of the question is the delineation of a genus. The 5% threshold of 16S has always been seen as a guide. I have never seen it rigidly enforced and it creates a conflict with phylogenetic theory as your example points out. Telling what is and ain't a genus ... you are right that a key is the tree structure.

1. **Tree structure and bootstrapping** ...

* A declared monophyly must be supported by a bootstrap \> 80%
* A declared paraphyly ... each paraphyletic branch must be supported by a bootstrap \> 80%

2. **16S versus other genes and tree structure** The tree structure needs confirming particularly as the bootstrap values 16S will often be low. The starting point is to consider the traditional MLST loci (I forget what they are but stuff like atpD, etc ...), build trees from them and :

* Don't combine 16S with protein genes (they evolve at different rates), although it could be done in a partitioned model
* MLST protein genes and core-house keeping genes can be combined without partitioning the data (there's a test, but it'll basically work).

Take the MLST genes repeat the tree and assess the bootstraps. Basically it will produce the same tree, but: a) confirms the tree structure; b) can provide stronger bootstrap values (which confirm the tree).

---

**5% stuff** All the above needs doing before any 5% stuff is considered. It is my personal view that using 5% 16S as a **mean value** is sort of bending the numbers, particularly in a paraphyly where the divergence is e.g. between 2 to 10%.

**Monophyly vs paraphyly** This where the tree structure is needed. Basically, a robust monophyly would be regarded as a taxonomic singular unit. In this case if its close to the 5% threshold genus rule:

* even if the pairwise distance between taxa within the monophyly is say 8% via 16S;
* and this is supported by other protein loci;
* and there's no suspicion this could represent a family.

In the above criteria I'd say its evidence of a genus. 10% is a bit big distance however.

If the tree structure is paraphyletic and the maximum different is say 10% between two taxa - one at the top of the paraphyly and the other as a most recent child ... then thats tricky. It could easily be two genera within the paraphyletic tree structure. For example, the vertebrates as paraphyletic, viz. amphibian, reptile, bird, edentates, marsupial, mammals and that represents 6 monophyletic Classes. Thus if that was my result I would leave open the possible interpretations: maybe this is a single genus; alternatively it could be two (or more) genera. The key to this is defining the nodes (internal branches) via bootstraps to be confident of the tree structure.

Generally a single robust paraphyletic tree structure, i.e. heavily imbalanced, where each paraphyletic branch is robustly supported would not be regarded as a single genus especially when 16S genetic differences is \>\> 5%. Taxonomically, this scenario would need further phylogenetic analysis (separate question) and classical bacterial taxonomists would want phenotypes before splitting a genus.

---

Reviewing the question against the answer - what I would do for the monophyly (diversity 10%) is split off the basal groups (at 10%) view robust bootstrap support. Then forward the in-group as the genus and the full monophyly as a putative genus. So that 10% is not the divergence unequivocally defining a single genus. So thats quite easy on the picture drawn. Thus the in-group is \<8% because 10% is a bit high.

It is harder to do this with a paraphyletic tree structure and thats a separate question on the advanced techniques needed to investigate this.

---
Tags: phylogenetics, statistics, 16rrna
---


# Question
Title: Calculating total number of Unique gene in every genome from ROARY pangenome output
I performed a ROARY from 100 genome. I got several file including gene\_presence\_absence.csv from the results. I want to know how many unique gene(unshared gene, based on blast homology, same gene (duplicate gene) may be unique if no homology is there between them) from presence absence csv file or other file. Kindly suggest me how to calculate total number of unique gene for each genome which is only present in one genome

# Answer
> 0 votes
Did you try query\_pan\_genome?

Ref. https://sanger-pathogens.github.io/Roary/

---
Tags: r, phylogeny, reference-genome
---


# Question
Title: How to pack up a Python script in the form a Dockerfile
I have a simple Python script with this header

```
from flask import Flask
app = Flask(__name__)

@app.route('/')
import os
import pandas as pd
from statsmodels.formula.api import ols
import statsmodels.api as sm

```

I have installed GitHub Desktop, Docker Desktop

I have created a folder in that I have

```
(.venv) (base) angel@Angels-MacBook-Pro docker % ls
Dockerfile      app.py          compose.yaml        requirements.txt

```

In `requirements.txt` I have

```
(.venv) (base) angel@Angels-MacBook-Pro docker % cat requirements.txt
blinker==1.6.2
click==8.1.6
colorama==0.4.6
Flask==2.3.2
itsdangerous==2.1.2
Jinja2==3.1.2
MarkupSafe==2.1.3
Werkzeug==2.3.6
gunicorn==21.2.0

```

This is the commands I have run so far

```
docker init

```

Then this one returns error

```
(.venv) (base) angel@Angels-MacBook-Pro docker % docker compose up --build
[+] Building 3.9s (14/14) FINISHED                                                                                docker:desktop-linux
 => [server internal] load .dockerignore                                                                                0.1s
 => => transferring context: 680B                                                                                0.1s
 => [server internal] load build definition from Dockerfile                                                                          0.1s
 => => transferring dockerfile: 1.63kB                                                                                0.1s
 => [server] resolve image config for docker.io/docker/dockerfile:1                                                                  2.3s
 => [server auth] docker/dockerfile:pull token for registry-1.docker.io                                                              0.0s
 => CACHED [server] docker-image://docker.io/docker/dockerfile:1@sha256:ac85f380a63b13dfcefa89046420e1781752bab202122f8f50032edf31b  0.0s
 => [server internal] load metadata for docker.io/library/python:3.9.12-slim                                                         1.3s
 => [server auth] library/python:pull token for registry-1.docker.io                                                                 0.0s
 => [server base 1/5] FROM docker.io/library/python:3.9.12-slim@sha256:364ada4cb4f943b603d330cbb3d26f689c8b54de51cb12bad1202de7dba3  0.0s
 => [server internal] load build context                                                                                0.0s
 => => transferring context: 1.39kB                                                                                0.0s
 => CACHED [server base 2/5] WORKDIR /app                                                                                0.0s
 => CACHED [server base 3/5] RUN adduser     --disabled-password     --gecos ""     --home "/nonexistent"     --shell "/sbin/nologi  0.0s
 => CACHED [server base 4/5] RUN --mount=type=cache,target=/root/.cache/pip     --mount=type=bind,source=requirements.txt,target=re  0.0s
 => [server base 5/5] COPY . .                                                                                0.0s
 => [server] exporting to image                                                                                0.0s
 => => exporting layers                                                                                0.0s
 => => writing image sha256:c207ecf929eb3dc4f172db22e56a19e7b89892fcdc31da80501cb3524e137f9d                                         0.0s
 => => naming to docker.io/library/docker-server                                                                                0.0s
[+] Running 1/1
 ✔ Container docker-server-1  Recreated                                                                                0.1s
Attaching to docker-server-1
docker-server-1  | [2023-09-22 11:09:20 +0000] [8] [INFO] Starting gunicorn 21.2.0
docker-server-1  | [2023-09-22 11:09:20 +0000] [8] [INFO] Listening at: http://0.0.0.0:8000 (8)
docker-server-1  | [2023-09-22 11:09:20 +0000] [8] [INFO] Using worker: sync
docker-server-1  | [2023-09-22 11:09:20 +0000] [9] [INFO] Booting worker with pid: 9
docker-server-1  | [2023-09-22 11:09:20 +0000] [9] [ERROR] Exception in worker process
docker-server-1  | Traceback (most recent call last):
docker-server-1  |   File "/usr/local/lib/python3.9/site-packages/gunicorn/arbiter.py", line 609, in spawn_worker
docker-server-1  |     worker.init_process()
docker-server-1  |   File "/usr/local/lib/python3.9/site-packages/gunicorn/workers/base.py", line 134, in init_process
docker-server-1  |     self.load_wsgi()
docker-server-1  |   File "/usr/local/lib/python3.9/site-packages/gunicorn/workers/base.py", line 146, in load_wsgi
docker-server-1  |     self.wsgi = self.app.wsgi()
docker-server-1  |   File "/usr/local/lib/python3.9/site-packages/gunicorn/app/base.py", line 67, in wsgi
docker-server-1  |     self.callable = self.load()
docker-server-1  |   File "/usr/local/lib/python3.9/site-packages/gunicorn/app/wsgiapp.py", line 58, in load
docker-server-1  |     return self.load_wsgiapp()
docker-server-1  |   File "/usr/local/lib/python3.9/site-packages/gunicorn/app/wsgiapp.py", line 48, in load_wsgiapp
docker-server-1  |     return util.import_app(self.app_uri)
docker-server-1  |   File "/usr/local/lib/python3.9/site-packages/gunicorn/util.py", line 371, in import_app
docker-server-1  |     mod = importlib.import_module(module)
docker-server-1  |   File "/usr/local/lib/python3.9/importlib/__init__.py", line 122, in import_module
docker-server-1  |     raise TypeError(msg.format(name))
docker-server-1  | TypeError: the 'package' argument is required to perform a relative import for '.venv.lib.python3.9.site-packages.eventlet.wsgi'
docker-server-1  | [2023-09-22 11:09:20 +0000] [9] [INFO] Worker exiting (pid: 9)
docker-server-1  | [2023-09-22 11:09:20 +0000] [8] [ERROR] Worker (pid:9) exited with code 3
docker-server-1  | [2023-09-22 11:09:20 +0000] [8] [ERROR] Shutting down: Master
docker-server-1  | [2023-09-22 11:09:20 +0000] [8] [ERROR] Reason: Worker failed to boot.
docker-server-1 exited with code 3

```

**yaml**

```
# Comments are provided throughout this file to help you get started.
# If you need more help, visit the Docker compose reference guide at
# https://docs.docker.com/compose/compose-file/

# Here the instructions define your application as a service called "server".
# This service is built from the Dockerfile in the current directory.
# You can add other services your application may depend on here, such as a
# database or a cache. For examples, see the Awesome Compose repository:
# https://github.com/docker/awesome-compose
services:
  server:
    build:
      context: .
    ports:
      - 8000:8000

# The commented out section below is an example of how to define a PostgreSQL
# database that your application can use. `depends_on` tells Docker Compose to
# start the database before your application. The `db-data` volume persists the
# database data between container restarts. The `db-password` secret is used
# to set the database password. You must create `db/password.txt` and add
# a password of your choosing to it before running `docker compose up`.
#     depends_on:
#       db:
#         condition: service_healthy
#   db:
#     image: postgres
#     restart: always
#     user: postgres
#     secrets:
#       - db-password
#     volumes:
#       - db-data:/var/lib/postgresql/data
#     environment:
#       - POSTGRES_DB=example
#       - POSTGRES_PASSWORD_FILE=/run/secrets/db-password
#     expose:
#       - 5432
#     healthcheck:
#       test: [ "CMD", "pg_isready" ]
#       interval: 10s
#       timeout: 5s
#       retries: 5
# volumes:
#   db-data:
# secrets:
#   db-password:
#     file: db/password.txt

```

**Dockerfile**

```
FROM python:3.9
WORKDIR /docker
COPY . /docker
RUN pip install -r requirements.txt
CMD ["python", "./app.py"]

```

# Answer
Not sure which tutorial you're following, but I think it's overly complicated. You do not need the compose file, and the Dockerfile has some issues:

First, specify what you want to copy into the container. That is `requirements.txt` and `app.py`. Then also, tell pip where the requirements file is. This works in terms of putting a container together. I am not testing the python script itself, this is on you.

```
$ cat Dockerfile
FROM python:3.9
WORKDIR /docker
COPY ["app.py", "requirements.txt" ,"/docker/"]
RUN pip install -r /docker/requirements.txt
CMD ["python", "./app.py"]

$ cat app.py
from flask import Flask
app = Flask(__name__)

@app.route('/')
import os
import pandas as pd
from statsmodels.formula.api import ols
import statsmodels.api as sm

$ cat requirements.txt
blinker==1.6.2
click==8.1.6
colorama==0.4.6
Flask==2.3.2
itsdangerous==2.1.2
Jinja2==3.1.2
MarkupSafe==2.1.3
Werkzeug==2.3.6
gunicorn==21.2.0

$ docker build -t atpoint/foo:1.0 .
[+] Building 5.8s (9/9) FINISHED
 => [internal] load build definition from Dockerfile                                                                                0.0s
 => => transferring dockerfile: 38B                                                                                0.0s
 => [internal] load .dockerignore                                                                                0.0s
 => => transferring context: 2B                                                                                0.0s
 => [internal] load metadata for docker.io/library/python:3.9                                                                                5.7s
 => [1/4] FROM docker.io/library/python:3.9@sha256:031b2a713079b2436e6c2102dbffbce109c9c2582ad0867c380e3c8d7fed477f                                                                                0.0s
 => [internal] load build context                                                                                0.0s
 => => transferring context: 64B                                                                                0.0s
 => CACHED [2/4] WORKDIR /docker                                                                                0.0s
 => CACHED [3/4] COPY [app.py, requirements.txt ,/docker/]                                                                                0.0s
 => CACHED [4/4] RUN pip install -r /docker/requirements.txt                                                                                0.0s
 => exporting to image                                                                                0.0s
 => => exporting layers                                                                                0.0s
 => => writing image sha256:6a6f1b08c8343608d2eacaa8463a1638491b504a3230ab00d45a7b32a41e1a43                                                                                0.0s
 => => naming to docker.io/atpoint/foo:1.0                                                                                0.0s

Use 'docker scan' to run Snyk tests against images to find vulnerabilities and learn how to fix them

```

Your container is now available as `atpoint/foo:1.0` which reads as username (atpoint), name of container (foo) and a version number / tag (1.0), all of which is free to choose. This is the most basic way of packing a container.
> 3 votes

---
Tags: python, docker, git
---


# Question
Title: Trimmomatic QC report shows drop in the reads and presence of overrepresented sequences
*This question was also asked on Biostars*

I am performing a *de novo* genome assembly using Illumina paired-end short reads, sequenced on a NovaSeq X by our collaborator at UCLA.

At present, I am in the stage of trimming the adapters. Here, you can have a look at the basic statistics and information on the adapter content obtained from the Fast QC report, for R1.

**FOR RAW READS**

I used Trimmomatic for trimming the adapter. The following is the **Trimmomatic Settings**

**ILLUMINACLIP:~/adapters/TruSeq3-PE.fa:2:30:10 MINLEN:36**

Here, you can see the basic statistics and adapter content of the Trimmed reads.

Here, the output was:

Both surviving: 566832403 Forward only surviving: 39244376 Reverse only surving: 0.00

Dropped reads: \<1%

Here you can go through the full FastQC report for:

Raw Reads

Trimmed Reads

**Now I have two questions:**

**Question 1**

Can I go ahead with the assembly process, for there is zero-adapter-presence in the reads? Should I mind the loss of reads?

**Question 2**

I see that there are over-represented sequences, both in read 1 and read 2. I doubt if I can leave them be, or I should trim them too. I would like to know if the over-represented sequences need to be tossed, or if they can be left untouched and proceed with the assembly. If I need to remove these sequences, can Trimmomatic do it?

If I were to trim the over-represented sequences, can I use an adapter "fasta" file as follows?

**The following are the over-represented sequences for R1.**

**The following are the over-represented sequences for R2.**

Can I have my fasta file as follows, to remove the over-represented sequences:

```
>PrefixPE/1
-------sequence1--------------
-------sequence2--------------
-------sequence3--------------
-------sequence4--------------
>PrefixPE/2
-------sequence1--------------

java -jar trimmomatic-0.30.jar PE s_1_1_sequence.txt.gz s_1_2_sequence.txt.gz
lane1_forward_paired.fq.gz lane1_forward_unpaired.fq.gz lane1_reverse_paired.fq.gz
lane1_reverse_unpaired.fq.gz ILLUMINACLIP:overRepresentedSeqs.fa:2:30:10  MINLEN:36

```

NB: Sorry about the lengthy post. Since the file sizes are large, it takes over a day for the Trimmomatic to finish the task, so that it is not feasible to go with multiple rounds of trial and error check.

# Answer
Given that it's a NovaSeq, these are almost certainly "no signal" reads as a result of the two-colour chemistry and should probably be excluded (or at least ignored).

A long polyG sequence is unlikely to affect mapping by much, because they're fairly rare in most genomes. Unfortunately you can't exclude them using Trimmomatic alone.

PolyG tail trimming can be done using fastp. I don't have any experience using fastp, so you'll need to read the documentation yourself.
> 1 votes

---
Tags: ngs, sequence-analysis, illumina
---


# Question
Title: 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?

# Answer
> 0 votes
Both of @haci's solutions are very cool:

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.

---
Tags: sequence-alignment, read-mapping, multiple-sequence-alignment, primer
---


title: For loop in Python goes wrong

text:
# Question
Title: For loop in Python goes wrong
I have mutation status (`0`,`1`) for ten genes in ten cancer 3D models

I have also gene expression in these 3D models, like

```
expr1 <- structure(list(GeneA_KO = c(0.0756383119639074, -0.162841965682723, -0.110703100204517, -0.457154616239399, -0.505815193885352, 0.366072775033448, 0.288650458935083, -0.0287646498640309, -0.690732824311905, 0.157366122620995, -1.29451546737684, 0.124167425460344, -0.168045313861751, 0.0186504943306465, -0.495880827400706, -0.206668627339351, -0.148609490335447, 0.066844529089804), GeneB_KO = c(-0.870157285953621, -0.35900872537626, -0.990406888128838, -0.120656981579292, -1.41427631791049, -0.866421095836781, -0.816936112146099, -0.894597236248507, -1.00084584207823, -0.850754488133456, -1.05314239761653, -1.07580675428372, -1.24708324665688, -1.04098460460899, -0.232388675985086, -1.48513703581865, -1.48373348870346, -1.76484058374354)), class = "data.frame", row.names = c("ModelM", "ModelO", "ModelS", "ModelA", "ModelI", "ModelC", "ModelT", "ModelH", "ModelE", "ModelR", "ModelA2", "ModelP", "ModelE2", "ModelU", "ModelT2", "ModelC2", "ModelS2", "ModelX"))
mut <- structure(list(ModelM = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelO = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelS = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelA = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelI = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelC = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelT = c(0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 0L), ModelH = c(1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 0L), ModelE = c(1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L), ModelR = c(1L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 0L), ModelA2 = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelP = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L), ModelE2 = c(1L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L), ModelU = c(1L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L), ModelT2 = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelI2 = c(1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelC2 = c(1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L),     ModelS2 = c(0L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 1L)), class = "data.frame", row.names = c("Gene1_mut", "Gene2_mut", "Gene3_mut", "Gene4_mut", "Gene5_mut", "Gene6_mut", "Gene7_mut", "Gene8_mut", "Gene9_mut", "Gene10_mut"))

```

`mut.csv` and `expr1.csv` exported for use in Python:

```
write.csv(mut, "mut.csv", quote = FALSE, row.names = FALSE)
write.csv(expr1[1:10,], "expr1.csv", quote = FALSE, row.names = FALSE)

```

I want to do a one-way ANOVA in Python so that calculate P-VALUE for each pair of mutation-expression in each pair of gene-model

This code returns error

```
import os
import pandas as pd
from statsmodels.formula.api import ols
import statsmodels.api as sm
m = pd.read_csv("mut.csv", sep=",")
m=pd.DataFrame(m)
m.info()
g = pd.read_csv("expr1.csv", sep=",")
g=pd.DataFrame(g)

lst = list(g.columns)
m_lst = list(m.columns)

dfs = pd.DataFrame({})
for Mut in m_lst:
    for Model in lst:
        df = pd.concat([m[[Mut]].reset_index(drop=True), g[[Model]].reset_index(drop=True)],axis=1)
        df = df.rename(columns={Mut: "mut", Model: "g"})
        model = ols('g ~ C(mut)', data=df).fit()
ValueError: endog has evaluated to an array with multiple columns that has shape (10, 18)

        aov_table = sm.stats.anova_lm(model, typ=2)
        aov_table[['Model']] = Model
        aov_table[['Mut']] = Mut
        dfs = dfs.append(aov_table)

```

The ideal output is something like

```
with pd.option_context('display.max_rows', None, 'display.max_columns', None):
    print(dfs)
            sum_sq   df         F    PR(>F)     Model      Mut
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO   ModelM
Residual  0.978877  8.0       NaN       NaN  GeneA_KO   ModelM
C(mut)    0.540949  1.0  7.375925  0.026417  GeneB_KO   ModelM
Residual  0.586718  8.0       NaN       NaN  GeneB_KO   ModelM
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO   ModelO
Residual  0.978877  8.0       NaN       NaN  GeneA_KO   ModelO

```

# Answer
> 1 votes
I can see what you're wanting to do, I can see bugs but the key bit of output is missing (below). It seems a bit weird to me the nested loop in context, but thats fine if you're happy.

The critical step is this

```
df = pd.concat([m[[Mut]].reset_index(drop=True), g[[Model]].reset_index(drop=True)],axis=1)

```

What I'm concerned about is the `pd.concat`. Are you sure `pd.concat` is the right command here? I thought `m` is one data table, e.g. mutation and `g` is another e.g. expression. `pd.concat` isn't going to work here, are you sure this isn't a `merge` or `join` command?

Anyway, the fitting is complaining the `.shape()` is wrong. Needs reshaping? Thats the error. Printing `df` at that point would be needed.

> ValueError: endog has evaluated to an array with multiple columns that has shape (10, 18)

To summarise normally the error you have received is about re-shaping the data. I've not put forward that solution here because,

* I am not sure what the data represents
* If I was doing this I wouldn't use `pd.concat`

It seems a lot of columns for a 1-way ANOVA, but could be right, I've not done 1-way ANOVA across 10 columns.

*Bug 1*

```
m = pd.read_csv("mut.csv", sep=",")
m=pd.DataFrame(m)
m.info()
g = pd.read_csv("expr1.csv", sep=",")
g=pd.DataFrame(g)

```

This should simply be:

```
m = pd.read_csv("mut.csv")
g = pd.read_csv("expr1.csv")

```

*Bug 2?*

```
dfs = pd.DataFrame({})
...
dfs = dfs.append(aov_table)

```

This will work given `aov_table` output is a pandas dataframe (I'm not at all sure it is). You need `aov_table.info()` to be sure. You can always place it in a dataframe however.

---
Tags: python, for-loop
---


# Question
Title: Would DESeq2 be appropriate or should I use another tool for differential gene expression?
Hay, this question is sort of an extension of previous post. We have sequenced extra cellular RNA from culture media. Since the number of cells were limited, we did not have enough material for biological or technical replicates, so each sample has only been sequenced once, and there are essentially three groups and each group has 2 subgroups, main groups are control, group A and group B. I think the table below should help -

| group | subgroup | Number of unique samples |
| --- | --- | --- |
| Control | Day 3 | ~30 |
| Control | Day 5 | ~5 |
| Treatment\_1 | Day 3 | ~30 |
| Treatment\_1 | Day 5 | ~50 |
| Treatment\_2 | Day 3 | ~40 |
| Treatment\_2 | Day 5 | ~20 |

We want to compare the gene expression (Or rather abundance of RNA) between samples based on the Day subgroup.

Since the same sample / culture batch has not been sequenced with replicates, and there is variance in number of samples per group, I was wondering if DESeq2 would be a valid option, as the documentation does mention that without replicates, the output from DESeq2 should be considered only for exploratory purpose.

If DESeq2 is not appropriate, which other tool would be better?

EDIT - I came across NOISeq which I will try once I get any output from DESEq2, if anyone has any other tools / papers which can potentially help me out, do let me know!

Thanks!

# Answer
> 1 votes
Definitely a fan of NOIseq :-) (long story) ... there's criticism it lacks power. What you are looking at is the significance of the observation against Type 1 error. The general criticism of non-parametrics is that they have low power, so few results. Personally, I'd use it and given how widely used DSeq2 is perform the same analysis via Dseq2. I'm not a fan of parametrics in this context, but in this field they dominate and DSeq2 is the best. Thus I'd use a parametric (DSeq2) and non-parametric NOIseq. The justification is you are using one of each. Personally, I do neither, I use machine learning to direct the thresholds of a parametric - but thats outside the public domain.

*The time series is more complicated.* Its essentially a 2-way problem: within a time point and between time points. Binary log difference is a good analysis between two time points in this situation.

I assume there was treatment between day 3 and day 5 for the treatment samples but not the control.

The issue of replicates. Personally what I would attempt to take a short cut with replicates: if Treatment 1 and 2 have a shared impact on the sample, then treating them as replicates and compare them against the control would be the first stage of analytics. This way if you go to the expense of producing the replicates described here you can be sure you've got a result before conducting further wet-lab work.

Clearly I don't know what Treatment 1 nor Treatment 2 are, but they could be drugs attacking the same metabolic pathway. Basically to do this you need a biological rationale for combining Treatments 1 and 2 to be considered replicates. I've attempted to re-describe this approach to biological analytics below the "-----"

The post-hoc analysis is more complex. The post-hoc analysis is a database operation assessing the relative binary log fold difference between time points for treatment versus control.

---

@haci likely represent the orthodox approach (e.g. at review level). What I would attempt is to provide biological justification that Treatment 1 and 2 were replicates. That depends on precise justification, i.e. this time series models the impact of generic interruption of this generic class of drug on pathway X, i.e. something that biochemically combines treatment 1 and 2.

The other analysis is a within time point analysis and if you did that between Treatment 1 and Treatment 2 for day 3 and day 5 separately (2 analyses), this would be used as justification if there ain't much different within a time point. **In addition** if treatment is post day 3, then Treatment 1, 2, and control are replicates, i.e. a triplicate (that's unequivocal). That provides a basis for moving forward.

However, The interaction @haci is describing is likely a formal model because this is a two-way problem. I rarely disagree but you've got to be careful here because its a time-series which aren't independent events and most "interaction" based models require independence. Thus day 5 is a result of day 3, therefore they are not independent. I've encountered this as an editor and it can go very wrong. Summary analysis between time series - those are independent, analysis between data points within a time-series are not independent and interaction mixes those up.

**Note** Time series experiments are always complex and performing wet-lab replicates would keep things simple (but obviously more expensive). So, if you've good analytics support or you're an analyst then you probably do it without further wet-lab, otherwise send it back to the wet-lab and then follow the recipe.

---
Tags: differential-expression, rna
---


# Question
Title: Align two different structures with respect to their active sites or heme cofactor
I have two enzymes, with heme cofactors. but their overall structure is very different.

I use Pymol regularly to visualize protein pdb structures.

I want to align both structures, but the structures don't align well due to very low homology. However, their active sites/substrate entrance channels may be more similar than the rest of the structure. so,

1. I want to align only the substrate channels from both enzymes and compare their structure. or
2. align both proteins only using the heme as a reference point. So both hemes from two enzymes align and then I can manually do the rest.

I am more interested in how the substrate would be accessible in both enzymes differently.

Thanks in advance.

Edit: Sorry for being a bit unclear. We can pick two proteins like cytochrome P450 BM3 with pdb id 6H1L or 6H1O. The second protein can be a similar peroxygenase UPO, e.g. 7pn4. Thanks.

# Answer
OK, trying to answer question 2, think I need to be sure that both enzymes have same identical heme cofactor and then find a rotation/translation matrix that applied to heme of protein A moved exactly on top of heme of protein B , than I have to move the entire protein backbone of protein A using that matrix, is it reasonable ? Am I getting it wrong ? With **PyMOL** I tried to use `Pair/Pair_fit` command with two heme proteins : `2d0t` (https://www.rcsb.org/structure/2d0t and `4zv8` (https://www.rcsb.org/structure/4ZV8) :

`pair hem_2d0t , hem_4zv8,` same as `pair_fit hem_2d0t , hem_4zv8`, where `hem_2d0t`,`hem_4zv8` are the heme of the two proteins \[`select model 2d0t and resn HEM` and `set_name sele, hem_2d0t` or `select hem_2d0t, model 2d0t and resn HEM`\].

See results, superposed heme:

Superposed Proteins :

As for question 1 I tried to answer using the tools available to me. I picked the two same proteins : `2d0t` (https://www.rcsb.org/structure/2d0t and `4zv8` (https://www.rcsb.org/structure/4ZV8). Using `parKVFinder` (https://github.com/LBC-LNBio/parKVFinder , https://www.softxjournal.com/article/S2352-7110(20)30319-8/fulltext) **pymol plugin** on both proteins stripped of all their HETATM, I calculated their cavities. As an example , here a png of `2d0t` cavities :

Then I extracted the residues surrounding the cavity using :

`select model 2d0t_merged near_to 5.0 of resn KAF` ; \[named `sele_2d0t`\]

`select model 4zv8_merged near_to 5.0 of resn KAB` ; \[named `sele_4zv8`\]

being `KAF` and `KAB` cavities described as non\_bounded atom in pdb format as output of `parKVFinder` .
I believe `parKVFinder` has is own output for the residues surrounding cavities, but wanted to try `PyMOL`

Example for `2d0t` :

Then I used

`cealign sele_4zv8 , sele_2d0t` (https://pymolwiki.org/index.php/Cealign) to align the residues of the two cavities. (`cealign` uses just CA of the selections provided).

Results , only cavities and surrounding residues :

the two heme protein cavities are not alignrd even if the HIS (2dt0, red) and CYS (4zv8, yellow) are on the same side :

NOTE:

2dt0 in magenta

4zv8 in cyan

I choose 2 different heme protein (could be more extreme though) , just because the question was open to broad attempts, could be interesting too see for example results on not too distant cytochromes
> 1 votes

---
Tags: pdb, pymol
---


title: Merge fragment reads

text:
# Question
Title: Merge fragment reads
I'm writing a function (python) that merges overlapping paired reads of a fragment, and a corresponding reference sequence. There's quite a lot of edge cases (e.g. when the reads disagree on indels) so I wonder if there's some standard library that does that?

Context - the input is pysam AlignedSegment DNA reads. The output is fragment, reference - both strings corresponding to the merged fragment and its corresponding reference.

e.g.

```
read1:           ACGTACGTACGT
read2:                CGTACGTACGT
expected output: ACGTACGTACGTACGT

read1:           ACGTANGTACGT
read2:               ACGTACGTACGT
expected output: ACGTACGTACGTACGT

```

this is just an illustration. The idea is to merge the overlapping DNA read sequences into the fragment and its corresponding reference using the information in the AlignedSegment (i.e. query\_sequence, reference\_start, query\_qualities, cigartuples, and get\_reference\_sequence())

Thanks.

# Answer
I'd suggest looking at similar tools that already exist, e.g. here is one that I hadn't heard of that lives on github. Just glancing at the code, it seems to use a scoring heuristic to decide how to derive a merge point, so it's using that to get around the absence of unambiguous points, probably at some rate of error.

BioPython has some functions that will create consensus sequences for you as well based on a pairwise alignment. (`Bio.Align.AlignInfo.SummaryInfo.dumb_consensus`, linked in that post.)

As a command line tool pandaseq does the same thing, I'd suggest looking at how they handle edge cases.

There are similar utilities in QIIME, NGMerge, BBMerge, \[illumina-utils\]- this is a pretty well-trodden path, I think all of those will be illuminating to look at in terms of reusing or writing original code.

I'm guessing it will be easier to merge FASTQs and then align, rather than align and then merge. (Note that these tools don't seem to operate on alignments, because that is trickier).

Alternately, if your real goal is to create a consensus genome sequence from your reference based on your alignments, there are a large number of tools that do that, notably `samtools` (see e.g. https://www.biostars.org/p/419927/). No new code needed, you just input your BAM file and reference sequence.
> 0 votes

---
Tags: pysam
---


# Question
Title: Generating CIGAR strings for MSA with reference
I have a list containing the results of a MSA:

```
msa = [
    '-G-ATG-CGGATACCG-', # reference sequence
    'CGTATCTCG--TACCGA',
    '-GTATCTCGT-GACCG-',
    'CG-AT--CGG-GACCGA',
    '---AT--CG--GACCGA',
    '-GTAT---GC-GACCG-',
    'CG-AT--CGG-GACCGA',
    'CGTAT--CG--GA CG-',
]

```

I would like to generate the CIGAR strings of those alignments (against the reference sequence); the expected result would be something like:

```
cigars = [
    '13M',
    '1I1M1I3M1I2M2D5M1I',
    '1M1I3M1I3M1D5M',
    '1I3M1D3M1D5M1I',
    '1D2M1D2M2D5M1I',
    '1M1I2M2D2M1D5M',
    '1I3M1D3M1D5M1I',
    '1I1M1I2M1D2M2D5M',
]

```

I don't know what would be a good approach for this. I managed to get the "aligned pairs" of each alignment but I don't think that'll help for determining the cigar strings:

```
# mapping = position in alignment => position in reference
mapping = {y: x for x,y in enumerate(y for y,nuc in enumerate(msa[0]) if nuc != '-')}

for aln in msa:
    aligned_pairs = [(y, mapping[z]) for y,z in enumerate(z for z,nuc in enumerate(aln) if nuc != '-') if z in mapping]
    print(aligned_pairs)

```

```
[(0, 0), (1, 1), (2, 2), (3, 3), (4, 4), (5, 5), (6, 6), (7, 7), (8, 8), (9, 9), (10, 10), (11, 11), (12, 12)]
[(1, 0), (3, 1), (4, 2), (5, 3), (7, 4), (8, 5), (9, 8), (10, 9), (11, 10), (12, 11), (13, 12)]
[(0, 0), (2, 1), (3, 2), (4, 3), (6, 4), (7, 5), (8, 6), (9, 8), (10, 9), (11, 10), (12, 11), (13, 12)]
[(1, 0), (2, 1), (3, 2), (4, 4), (5, 5), (6, 6), (7, 8), (8, 9), (9, 10), (10, 11), (11, 12)]
[(0, 1), (1, 2), (2, 4), (3, 5), (4, 8), (5, 9), (6, 10), (7, 11), (8, 12)]
[(0, 0), (2, 1), (3, 2), (4, 5), (5, 6), (6, 8), (7, 9), (8, 10), (9, 11), (10, 12)]
[(1, 0), (2, 1), (3, 2), (4, 4), (5, 5), (6, 6), (7, 8), (8, 9), (9, 10), (10, 11), (11, 12)]
[(1, 0), (3, 1), (4, 2), (5, 4), (6, 5), (7, 8), (8, 9), (9, 10), (10, 11), (11, 12)]

```

How can I generate the CIGAR strings efficiently?

# Answer
If you loop through the characters of two aligned sequences **in parallel** then you can easily determine the operations that compose the CIGAR. Then you just have to "collapse" the generated list of operations:

```
import itertools

ref = msa[0]

for aln in msa:
    cig = []

    for i,nuc in enumerate(aln):
        if nuc != '-':
            if ref[i] != '-':
                cig.append('M')
            else:
                cig.append('I')
        elif ref[i] != '-':
                cig.append('D')

    cig = [(len(list(it)), op) for op,it in itertools.groupby(cig)]

    print( ''.join(str(n) + op for n,op in cig) )

```

```
13M
1I1M1I3M1I2M2D5M1I
1M1I3M1I3M1D5M
1I3M1D3M1D5M1I
1D2M1D2M2D5M1I
1M1I2M2D2M1D5M
1I3M1D3M1D5M1I
1I1M1I2M1D2M2D5M

```
> 4 votes

---
Tags: python, cigar
---


# Question
Title: Kraken2 Standard Database failing to build (unexpected FTP path)
I am attempting to build just the standard Kraken2 Database by using the following command:

```
kraken2-build --standard --threads 24 --db $DBNAME

```

I am returned the following error a few minutes building it:

```
rsync_from_ncbi.pl: unexpected FTP path (new server?) for
https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/762/265/GCF_000762265.1_ASM76226v1

```

Anyone know why this is happening?

# Answer
I have experienced the same issue yesterday and after trying the `--use-ftp` argument without luck (this is proposed as an alternative in the forums), I have downloaded the database from Kraken 2, KrakenUniq and Bracken indexes (ref).
> 1 votes

---
Tags: ngs, metagenome, ncbi
---


# Question
Title: How can I develop a snakemake rule for all files in multiple directories?
I am a beginner in snakemake and I'm trying to develop a working flow.

Imagine I have multiple folders and in each folder I have a tumour and normal sample files:

```
SRC120/
      SRC120-T1.fastq
      SRC120-N1.fastq

SRC121/
      SRC121-T1.fastq
      SRC121-N1.fastq

```

and so on (in total I have more than 100 SRC folders).

I am trying to make a snakemake rule to make a symlink of these fastq files.

What I have written so far is:

```
import os

SAMPLES=os.listdir("/analysis/Anna/snakemake-demo/data/samples_fastq-test1/")

def getFastqFilesForTissu(wildcards):
    fastqs = list()
    # Beware no other files than fastqs should be there
    for s in os.listdir("/analysis/Anna/snakemake-demo/data/samples_fastq-test1/"+wildcards.sample):
        fastqs.append(os.path.join("/analysis/Anna/snakemake-demo/data/samples_fastq-test1/",wildcards.sample,s))
    return fastqs

rule all:
    input:
        expand("results/{sample}/{sample}.fastq", sample=SAMPLES)

rule foo:
    input:
        fastq_files = getFastqFilesForTissu
    output:
        "results/{sample}/{sample}.fastq"
    shell:
        'ln -s {input} {output}'

```

This works if I have one file in each folder and when the file name matches with the folder's name (i.e., when there is no T1 or N1 suffix).

Is there any way I could modify my rule to get it working for all samples in SRC directories? Or a better way to handle this issue?

# Answer
I suggest using package glob, which supports unix-style pattern search. For instance finding all .py files in my repo:

```
import glob
all_python_files = glob.glob(r'c:\personal\coderepo\**\*.py', recursive=True)

```

Which outputs full paths:

> 'c:\personal\coderepo\SimpleHTR\data\checkDirs.py',
> 'c:\personal\coderepo\\....

I typically like to input the glob.glob to Pandas dataframes:

```
all_files = glob.glob(f'{base_folder}/*')
all_files = pd.DataFrame(all_files, columns=['full_path'])
all_files['my_filename'] = all_files['full_path'].apply(lambda x: os.path.basename(x))

#do whatever filtering is required:
all_files = all_files[all_files['my_file_name'].str.contains('regex_here')]

```

Then perform what ever I need to do in loop:

```
for idx, _ in all_files.iterrows():
    do_something(_['full_path'])

```
> 2 votes

---
Tags: python, snakemake
---


title: Parameter optimisation of minimap2

text:
# Question
Title: Parameter optimisation of minimap2
I'm leveraging **minimap2** to map select genes from short-read microbial fastq metagenomics zip files.

```
minimap2 --heap-sort=yes --splice-flank=no metagenome.fastq.gz fastaDB.fas

```

What I noticed is the low tolerance of mismatch when I need to maximise mismatch because microbes have loads of genetic diversity. Thus the collection of references sequences are never really a reference, in fact the "microbe reference genome" is a misnomer, because *in situ* diversity of microbes is huge.

**Question** How do you relax the mapping stringency criteria for minimap2? If you are familiar with Blast parameters this would help alot here. I know blast but I've never parameterised it.

The options are this are:

*Alignment options*

> -A Matching score \[2\]

> -B INT Mismatching penalty \[4\]

> -O INT1\[,INT2\] Gap open penalty \[4,24\].

> -s INT Minimal peak DP alignment score to output \[40\].

> --score-N INT Score of a mismatch involving ambiguous bases

Anything else that is relevant.

Note, -f is a famous option, not sure it helps here.

I'm going to code up a hyperparameterisation unit-test, but understanding the parameters first would really, really help.

Alternative approaches welcome as well.

---

*Plugging minimap2*. I code high performance (hadoop-style, back-tracking, graph, traverse tall buildings in a stride etc ...). The coder that wrote minimap has knowledge beyond, i.e. very fast run times directly from gzip file. I assume they are indexing then chunking the raw zip file then unzipping each 'chunk' JIT (just in time).

---

**Thanks @MaximilianPress great answer**. I'll fill in the backstory.

I'm trapping total genetic diversity across Genbank and the metagenomes. My protocol to do both:

* Blasting out Genbank and then
* piping the output across metagenomes via in this case minimap2.

Hence speed is key. So I ain't replacing Blast with minimap2 but leveraging minimap2 via the output of blast.

*However* on reflection of @MaximilianPress's discussion the answer is to generate a gold standard of diversity to assess alternative approaches.

# Answer
> 2 votes
## First of all, use short read option

The short read preset (`-x sr`) for minimap2 is heavily recommended for doing such mapping:

`minimap2 -x sr`

## Heng Li on minimap2 for related genomes

Heng Li (minimap2 author) has some suggestions on a somewhat related github issue (though it is now 4 years old):

1. Minimap2 is not really intended for this purpose. For short reads he recommends stampy, but it is clearly quite old (python2.7) and does not seem to have much recent developer history.
2. Failing that, he suggests an option set such as follows:

> If you want to give minimap2 a try anyway, you can use the following setting: `minimap2 -axsr -A2 -B4 -O4,24 -E2,1` If minimap2 is not sensitive enough, reduce -k and -w at the cost of mapping speed. Reducing -w should have relatively lower impact on speed, but reducing -k should be more sensitive.

I would suggest engaging with Heng, he is quite responsive in general and may have more suggestions.

## Other aligners for metagenomics

I would suggest that you look for more possibilities for alignment, or possibly not using alignment in the first place. However, this is very sensitive to your ultimate application, which I'm not sure that I understand from your post.

I was pleased to see that others have reported good results using minimap2 for metagenomics, which I was skeptical about, having been out of the field for a bit. Usually in metagenomics there is an expected tradeoff between sensitivity and specificity for the reasons that you state (this paper is old but it uses BLAST, so it is a sort of best case for sensitivity on unrelated genomes at the level of divergence that they test).

There are of course new shotgun metagenomic aligners coming out all the time, here is one recent example that I know nothing about but which may be interesting as a comparison.

If you are really looking for a rich parameter set that allows low-stringency mappings, I would suggest going back to BLAST, though it will be much slower. minimap2 and related algorithms are to some extent dependent on the ability to find short seed matches, so there will necessarily be less than perfect sensitivity in your application.

# Answer
> 1 votes
Just to follow-up on the performance of **minimap2**: the solution @MaximillianPress described works. Thousands of hits are obtained for comparatively obscure genes (some will be duplicates). Parameterisation gives a 2-log10 increase (astonishing).

The issue is that its:

* relatively slow - precisely as described by @MP;
* very, very sensitive to parameterisation, just a single integer changes the output by a couple of thousand hits.

I can see why there was concern about using it - such notable sensitive is an issue - and alternatively its clear why other investigators have had great success with it - hit the "sweet-spot" and its a bonanza. For my purposes it's fine - it's not a quantitation of complete genetic diversity, as long as I capture a reasonable proportion of diversity that is fine. Clearly, they'll be times the "sweet-spot" is missed but thats still okay (famous last words :-) ).

---
Tags: fastq, blast, read-mapping, minimap2, parameters
---


# Question
Title: ValueError: count not convert string to float
I'm trying to use a Python script to plot with a rolling window using pandas and seaborn. This code worked for the longest time but now it's giving me an error that I don't know how to fix.

Here is the code I'm using:

```
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

heterozygosity_df = pd.read_csv("file.tsv", sep="\t", header=None).set_axis(['chr', 'pos', 'het'], axis=1, copy=False)

heterozygosity_df['het'] = heterozygosity_df['het'].transform(lambda x: -1*(x-1))

windows = [25, 50, 100, 200]
steps = [5]

for window in windows:
    for step in steps:
        heterozygosity_ratio = {}
        for chr_heterozygosity_df in heterozygosity_df.groupby('chr'):
            heterozygosity_ratio[chr_heterozygosity_df[0]] = chr_heterozygosity_df[1].set_index('pos').rolling(window=window, step=step, center=True).mean()
...

```

There is more after this but that's just to configure the resulting plot with matplotlib. The error I get when I use the above code is this:

```
Traceback (most recent call last):
  File "/opt/homebrew/lib/python3.11/site-packages/pandas/core/window/rolling.py", line 368, in _prep_values
    values = ensure_float64(values)
             ^^^^^^^^^^^^^^^^^^^^^^
  File "pandas/_libs/algos_common_helper.pxi", line 42, in pandas._libs.algos.ensure_float64
ValueError: could not convert string to float: 'chrI'

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

Traceback (most recent call last):
  File "/opt/homebrew/lib/python3.11/site-packages/pandas/core/window/rolling.py", line 485, in _apply_blockwise
    arr = self._prep_values(arr)
          ^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/homebrew/lib/python3.11/site-packages/pandas/core/window/rolling.py", line 370, in _prep_values
    raise TypeError(f"cannot handle this type -> {values.dtype}") from err
TypeError: cannot handle this type -> object

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

Traceback (most recent call last):
  File "plot_heterozygosity_window.py", line 29, in <module>
    heterozygosity_ratio[chr_heterozygosity_df[0]] = chr_heterozygosity_df[1].set_index('pos').rolling(window=window, step=step, center=True).mean()
                                                     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/homebrew/lib/python3.11/site-packages/pandas/core/window/rolling.py", line 2111, in mean
    return super().mean(
           ^^^^^^^^^^^^^
  File "/opt/homebrew/lib/python3.11/site-packages/pandas/core/window/rolling.py", line 1507, in mean
    return self._apply(window_func, name="mean", numeric_only=numeric_only)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/homebrew/lib/python3.11/site-packages/pandas/core/window/rolling.py", line 615, in _apply
    return self._apply_blockwise(homogeneous_func, name, numeric_only)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/homebrew/lib/python3.11/site-packages/pandas/core/window/rolling.py", line 487, in _apply_blockwise
    raise DataError(
pandas.errors.DataError: Cannot aggregate non-numeric type: object

```

I have tried converting the 'chr' column directly from a string to a float: `heterozygosity_df['chr'] = heterozygosity_df['chr'].astype(float)`

I have tried converting the entire dataframe to a float: `heterozygosity_df = heterozygosity_df.astype(float, errors = 'raise')`

But everything keeps giving me the `ValueError: could not convert string to float: 'chrI'`

Does someone know an easy way to fix this? Like I said, I have been using this code for months, and it's always worked - have they deprecated something in the pandas package to make it so you can't read in columns with strings and integers?

**EDIT**:

For reference the dataframe is structured like this:

```
chrI  123  1
chrI  234  0
chrI  237  1
chrII 456  1
chrII 567  1
chrII 678  1
...
chrXVI 123 1

```

Where the first column designates the chromosome, the second column is the position, and the third column is the heterozygosity indicator.

# Answer
The following would be useful,

```
print (heterozygosity_df.columns)

```

This looks like static typing issue within pandas. What I suspect is `heterozygosity_df.['chrI']` is a column in the dataframe. What I think has happened is there's a mix of strings and floats within this column. pandas has set this as a "string" but you are wanting to perform numerical operations. Thus the solution is simply

```
print(heterozygosity_df.dtypes) # this should state 'chrI' is a "category" or "string"
heterozygosity_df['chrI'] = heterozygosity_df['chrI'].astype(float)
print(heterozygosity_df.dtypes)

```

If you have multiple changes the syntax is

```
heterozygosity_df = heterozygosity_df.astype({'chrI':'float', 'egColumn':'category'})

```

I suspect there will be other errors, e.g. the header is the first row of the data column. This is because again for pandas to automatically assign a column to a "string" means there must be a string value within the column.

---

From the comments. I see whats happening. The easy solution is ...

```
replacement = {
    "chrI": 1,
    "chrII": 2,
    "chrIII": 3,
    "chrIV": 4,
    ....
}

heterozygosity_df['chr'] = heterozygosity_df['chr'].str.replace(replacement, regex=True)

```

---

From the comments ... good it works ... this is what I would have personally done ..

```
replacement = {
    "chrI": 1,
    "chrII": 2,
    "chrIII": 3,
    "chrIV": 4, # continue for all chromosomes
}

heterozygosity_df = pd.read_csv("file.tsv", sep="\t", header=None).set_axis(['chr', 'pos', 'het'], axis=1, copy=False)
heterozygosity_df['chr'] = heterozygosity_df['chr'].str.replace(replacement, regex=True).astype('int')

```

Normally `str` should be in place because when `'chr'` is imported its an object. However, if it works thats the only thing that counts.
> 1 votes

# Answer
Ah! I finally got it to work! I used your command but I had to change it up a little bit to work with my data. The final command I used is this:

```
replacement = {
    "chrI": 1,
    "chrII": 2,
    "chrIII": 3,
    "chrIV": 4,
    "chrV": 5,
    "chrVI": 6,
    "chrVII": 7,
    "chrVIII": 8,
    "chrIX": 9,
    "chrX": 10,
    "chrXI": 11,
    "chrXII": 12,
    "chrXIII": 13,
    "chrXIV": 14,
    "chrXV": 15,
    "chrXVI": 16,
    "chrmt": 17
}

heterozygosity_df['chr'] = heterozygosity_df['chr'].replace(replacement, regex=False)

```

And with that I was able to generate plots showing all of the chromosomes! (I have to figure out how to change the axis labels from 1,2,3,4...10 back to the chromosome but that's a future problem). Thank you so much for all your help!!
> 2 votes

---
Tags: python, errors, pandas, seaborn
---


title: tensorQTL interaction issue

text:
# Question
Title: tensorQTL interaction issue
I hope that the question has its place on this board.

I am trying to run the mapping of cis eQTL using an interaction term with tensorQTL. I am running into the following index error when running the software:

```
>>> cis.map_nominal(genotype_df, variant_df,
...                 phenotype_df,
...                 phenotype_pos_df,
...                 prefix = prefix,
...                 covariates_df=covariates_df,
...                 interaction_df=interaction_s,
...                 maf_threshold_interaction=0.05,
...                 run_eigenmt=True,
...                 output_dir=output,
...                 write_top=True,
...                 write_stats=True)
cis-QTL mapping: nominal associations for all variant-phenotype pairs
  * 569 samples
  * 24665 phenotypes
  * 68 covariates
  * 69763935 variants
  * including 1 interaction term(s)
    * using 0.05 MAF threshold
  * cis-window: ±1,000,000
  * checking phenotypes: 24665/24665
  * Computing associations
    Mapping chromosome 1
    processing phenotype 137/24665Traceback (most recent call last):
  File "/QRISdata/Q5776/software/tensorqtl/tensorqtl/core.py", line 197, in calculate_interaction_nominal
    Xinv = torch.matmul(torch.transpose(X_t, 1, 2), X_t).inverse() # ng x (1+2*ni) x (1+2*ni)
torch._C._LinAlgError: linalg.inv: (Batch element 2704): The diagonal element 3 is zero, the inversion could not be completed because the input matrix is singular.

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/QRISdata/Q5776/software/tensorqtl/tensorqtl/cis.py", line 294, in map_nominal
    res = calculate_interaction_nominal(genotypes_t, phenotype_t.unsqueeze(0), interaction_t,
  File "/QRISdata/Q5776/software/tensorqtl/tensorqtl/core.py", line 200, in calculate_interaction_nominal
    i = int(re.findall('For batch (\d+)', str(e))[0])
IndexError: list index out of range

```

However, when I run the same command but only on 50 genes, the code runs properly.

```
>>> cis.map_nominal(genotype_df, variant_df,
...                 phenotype_df.iloc[1:50],
...                 phenotype_pos_df.iloc[1:50],
...                 prefix = prefix,
...                 covariates_df=covariates_df,
...                 interaction_df=interaction_s,
...                 maf_threshold_interaction=0.05,
...                 run_eigenmt=True,
...                 output_dir=output,
...                 write_top=True,
...                 write_stats=True)
cis-QTL mapping: nominal associations for all variant-phenotype pairs
  * 569 samples
  * 49 phenotypes
  * 68 covariates
  * 69763935 variants
  * including 1 interaction term(s)
    * using 0.05 MAF threshold
  * cis-window: ±1,000,000
  * checking phenotypes: 49/49
  * Computing associations
    Mapping chromosome 1
    processing phenotype 49/49
    time elapsed: 0.44 min
    * writing output
done.

```

The interaction factor is a categorical variable (0/1), I am not sure why the issue is being raised and how to fix

I would appreciate any help/direction to go forward.

Thanks

# Answer
> 1 votes
Writing a solution for anyone interested. The issue stems from a dosage error, some variant contains only heterozygous individuals and no homozygous ones needed for QTL calculation.

I used plink2 to calculate genotypes for all variant using the following command:

```
plink2  --pfile PATH \
     --geno-counts  \
     --out OUT

```

I then removed variant that contains less than 10 Homozygous invidivuals (HOM\_REF\_CT \< 10) using plink then was able to run tensorQTL interaction mode.

Hope it helps.

---
Tags: genomics, genetics, qtls
---


title: Nanopore Flongle vs MinION

text:
# Question
Title: Nanopore Flongle vs MinION
How do the Oxford Nanopore Flongle flow cells perform compared to the standard MinION flow cells, when it comes to reliability and error rates? I found this comparison article from last year, from where the below graph is taken, but they used the old chemistry. How is the situation today with the new R10.4.1 nanopores and Kit 14 chemistry? Does anyone have first-hand experience, or a good reference?

My use case would mainly be metabarcoding, i.e. amplicon sequencing. If anyone has experience of Flongle flow cells for specifically metabarcoding, I would be happy to hear your thoughts.

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

Here are the yield assumptions I use for flow cells, which are basically the lower ranges that we've seen when using them:

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

For Flongle, this has been fine for us to do sequencing and assembly of 12 plasmids on one flow cell using the rapid barcoding kit, and the improved accuracy from the Kit14 + R10.4.1 combination makes reliable assembly easier (especially where homopolymers are involved). I expect it will work similarly well for amplicon sequencing.

Error rates are pretty similar across platforms, because it's the pore that matters, not the package. Here's a comparison plot of mapped accuracy vs predicted accuracy (from the basecaller) that I did recently on a TrackIt 1kb Plus DNA ladder, with bead cleanup using LFB, run on a PromethION flow cell with the V14 ligation kit (sequenced on a P2 Solo), then called using the new bacterial model. The reads represented in this plot have been filtered to only keep (at most) the longest 2000 matches for each target. The run produced 1,128,671 reads (1.95 Gb) over 1.5 hours:

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

The predicted accuracy is taken from the basecaller (i.e. fastq quality string), and is calculated by summing up the quality values after converting to a linear representation, then dividing by the number of bases and converting back to a log representation:

```
sub getMeanQual {
  my ($qual) = @_;
  if(length($qual) == 0){
    return(0);
  }
  my $qBase = 33;
  my @qspl = split("", $qual);
  my $qualTotal = 0;
  my $qualCount = 0;
  foreach my $qual (@qspl) {
    $qualTotal += 10 ** (-(ord($qual) - $qBase) / 10);
    $qualCount++;
  }
  return(int(-10 * log($qualTotal / $qualCount) / log(10)));
}

```

According to the trained models that LAST generated, there's not much different between the called reads between Flongle and PromethION flow cells. Here's the PromethION base mismatch matrix (q45 SNPs):

```
# substitution percent identity: 99.9971
#last -Q 1
#last -t4.38507
#last -a 23
#last -A 19
#last -b 2
#last -B 4
#last -S 1
# score matrix (query letters = columns, reference letters = rows):
       A      C      G      T
A      6    -45    -35    -46
C    -45      6    -45    -44
G    -33    -45      6    -46
T    -46    -45    -46      6

```

And here's the Flongle base mismatch matrix (q41 SNPs):

```
# substitution percent identity: 99.9922
#last -Q 1
#last -t4.36483
#last -a 20
#last -A 20
#last -b 4
#last -B 5
#last -S 1
# score matrix (query letters = columns, reference letters = rows):
       A      C      G      T
A      6    -45    -36    -46
C    -45      6    -44    -45
G    -32    -45      6    -46
T    -46    -45    -46      6

```

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

---
Tags: nanopore, minion, barcode
---