# Question
Title: What is the best way to compare differences in DEG expression between two populations?
I am interested in differential response to a drug in the "RED" group versus the "GRN" group. I mean, I need to compare differences in DEG expression between two populations (RED vs GRN)

I have 3 pairs of control vs treated samples. 3 pairs for each population.

Here is an example of results for a gene named "gen1".

```
Gname   RED_ctl1  RED_drg1  RED_ctl2  RED_drg2  RED_ctl3  RED_drg3
gen1           4        10         5        10         6         9
Gname   GRN_ctl1  GRN_drg1  GRN_ctl2  GRN_drg2  GRN_ctl3  GRN_drg3
gen1           4         9         5        10         4         9

```

In the end, I expect to plot a boxplot that would show me (or not) the difference between RED and GRN groups in terms of DEG expression. I expect to find a group (either RED or GRN) with DEGs that respond to drugs much more (or less) in terms of expression.

Since this is my first time doing this analysis, I would like to know what is the classical way of analyzing that kind of data? What are the steps?

A related question was presented here on Bioinformatics SE

# Answer
A t-test or eyeballing a box plot might be okay for something very informal, but that's not the statistically right way to do it.

If this is based on RNASeq data, the answer is clear. You need DESeq2 or EdgeR. Read their tutorials, and pick one to use.

If this is microarray, limma is probably best.

If it's something else, you'll need to read some papers to see what people use.
> 2 votes

# Answer
It doesn't matter what data you are dealing with classical statistics can handle this data and non-parametric are very versatile, but there is one notable limitation.

If this is strictly paired data (which what I suspect):

```
GRN 1 | GRN 2 | GRN 3
RED 1 | RED 2 | RED 3

```

This would be a two-way ANOVA and is a very powerful test it is fulfils the criteria (below). Interpreting the output needs a bit of background reading.

However, if it is simply GRN vs RED then it is one-way ANOVA (or T-test) and Kruskal-Wallis test (non-parametric). In addition, if you are comparing

```
RED1
RED2
RED3
GRN1
GRN2
GRN3

```

Thus 6 data sets a T-test doesn't work and its one-way ANOVA or KW test.

*ANOVA* The Fmax test is needed for ANOVA analysis to assess the homogeneity of variance.

Additional ANOVA criteria

* The test works best when the number of samples drawn from each population is the same (check in this case). This is not an essential requirement but if the sampling between groups is heavily that's not good.
* The underlying populations are normally distributed (same is true of the T-test). If Fmax says 'yes' skip this bit.

I honestly don't know if Fmax is available in Python, it is definitely not `numpy.fmax`. It is available in R:

```
hartleyTest(x, ...)

```

where *x* is the vector. You can do it manually.

For 1-way ANOVA

```
from scipy.stats import f_oneway

f_oneway(GRN123, RED123)
f_oneway(GRN1,GRN2,GRN3,RED1,RED2,RED3) # This makes it like the second example

```

The data is within a list e.g. `GRN1 = [1,2,3,1,1,14,1]`

Easy.

Two-way ANOVA is complicated

```
import statsmodels.api as sm
from statsmodels.formula.api import old

model = ols('treatment ~ C(RED1) + C(GRN1) ... : C(GRN1)', data=df).fit() # see note
sm.stats.anova_lm(model, typ=2)

```

You should be able to circumvent the ols - because its complicated. It is however useful, particularly is Fmax fails.

Here for details.

You might substitute your one-way ANOVA for a T-test in this instance, but it depends how you structure you comparison.

*If Fmax fails* and you were wanting to perform a one-way ANOVA you shift to the non-parametric Kruskal-Wallis test. In python,

```
from scipy import stats

stats.kruskal(GRN123, RED123)

```

Work around if its paired data ...

```
stats.kruskal (GRN1, RED1)
stats.kruskal (GRN2, RED2)
....

```

Its a bit cheeky but you can cover yourself with Bonferoni correction and strictly Wilcoxan is better here.

---

So its easy.

**HOWEVER** if your data is paired and fails an Fmax test, strictly you go to the Friedman test, this is the non-parametric version of the 2-way ANOVA (by Milton Friedman). Now for some reason no-one ever does this, but I don't remember what the reason why. Thus your options are:

* in that situation you could just use the ols regression model fitted into the ANOVA - but you need to think hard about his model because it needs to present nicely paired data to ANOVA.
* You could of course do the Friedman test and its up to the reviewer to understand the theoretical criticism of this (I don't remember it and generally the quantitive standards in molecular peer-review isn't great)
* You could use multiple KW

**Alternative** the 'correct' way to get ANOVA working when Fmax says 'no' is transformation. This gets complicated however and is one step before a custom model.

---

Final points

1. **Exception** If this is microarray data - there is an issue about standardisation and that is complicated (comparing different microarray slides against experimental variation).
2. What algorithms such as DSeq2 do will likely do is apply a fixed transformation to the data, but I don't know about standardisation. **Non-parametric statistics will circumvent the requirement for transformation** ... their whole point is you don't transform and must not transform - must be raw data.
3. Again non-parametrics become difficult in replicating 2-way ANOVA ... there are several approaches here and of course there is nothing to stop you using it and leave to someone else explain the limitation.

**Note** the Wilcoxon is the non-parametric equivalent of the T-test, which I think we covered in a previous post (these are mice right?).

---

What would I do? Personally, if a non-parametric fits I use it every time ... its simple, everyone understand it and a significant result is not in dispute. The reality is you are probably into a skewed distribution. Furthermore, no-one gets worried about a few KW tests (type 2 error) and you can always deploy bonferoni's correction.

I should understand the limitation of Friedman test and I might be working under a false assumption. If a simple transformation results in good Fmax I'd use ANOVA (if raw data isn't homogeneous). Beyond its statistical modelling and that's complicated.
> 1 votes

---
Tags: differential-expression
---


title: Problems when using vcf2maf

text:
# Question
Title: Problems when using vcf2maf
I am trying to use a commonly recommended package named vcf2mafto convert vcf file to maf file and is reliant on Ensembl's Variant Effect Predictor (VEP) here. However, if I follow the guidelines shown on the github, I would get:

```
perl vcf2maf.pl --ref-fasta hg38.fa --input-vcf chr21.vcf --output-maf chr21.maf
dyld[33981]: Library not loaded: @rpath/libcrypto.1.0.0.dylib
  Referenced from: /usr/local/anaconda3/bin/samtools
  Reason: tried: '/usr/local/anaconda3/bin/../lib/libcrypto.1.0.0.dylib' (no such file), '/usr/local/anaconda3/bin/../lib/libcrypto.1.0.0.dylib' (no such file), '/usr/local/lib/libcrypto.1.0.0.dylib' (no such file), '/usr/lib/libcrypto.1.0.0.dylib' (no such file)
dyld[33982]: Library not loaded: @rpath/libcrypto.1.0.0.dylib
  Referenced from: /usr/local/anaconda3/bin/samtools
  Reason: tried: '/usr/local/anaconda3/bin/../lib/libcrypto.1.0.0.dylib' (no such file), '/usr/local/anaconda3/bin/../lib/libcrypto.1.0.0.dylib' (no such file), '/usr/local/lib/libcrypto.1.0.0.dylib' (no such file), '/usr/lib/libcrypto.1.0.0.dylib' (no such file)
dyld[33983]: Library not loaded: @rpath/libcrypto.1.0.0.dylib
  Referenced from: /usr/local/anaconda3/bin/samtools
  Reason: tried: '/usr/local/anaconda3/bin/../lib/libcrypto.1.0.0.dylib' (no such file), '/usr/local/anaconda3/bin/../lib/libcrypto.1.0.0.dylib' (no such file), '/usr/local/lib/libcrypto.1.0.0.dylib' (no such file), '/usr/lib/libcrypto.1.0.0.dylib' (no such file)
dyld[33984]: Library not loaded: @rpath/libcrypto.1.0.0.dylib
  Referenced from: /usr/local/anaconda3/bin/samtools
  Reason: tried: '/usr/local/anaconda3/bin/../lib/libcrypto.1.0.0.dylib' (no such file), '/usr/local/anaconda3/bin/../lib/libcrypto.1.0.0.dylib' (no such file), '/usr/local/lib/libcrypto.1.0.0.dylib' (no such file), '/usr/lib/libcrypto.1.0.0.dylib' (no such file)
ERROR: You're either using an outdated samtools, or --ref-fasta is not the same genome build as your --input-vcf. at vcf2maf.pl line 410.

```

I think the reference is the correct one, and I guess the error might be because of the VEP installation (as I came across a lot of bugs when following their gist).

As a result, I tried to use the Conda installation and use the codes:

```
conda install -c bioconda vcf2maf

```

However, after the package was successfully installed, I could not call it in the command line.

```
(base) user-MacBook-Pro:~ scottxu$ conda install -c bioconda vcf2maf
Collecting package metadata (current_repodata.json): done
Solving environment: done

# All requested packages already installed.

(base) user-MacBook-Pro:~ user$ vcf2maf
-bash: vcf2maf: command not found

```

Or, if I want to use the Perl script: I will get:

```
(base) user-MacBook-Pro:~ user$ perl vcf2maf.pl
Can't open Perl script "vcf2maf.pl": No such file or directory

```

Could anyone give me suggestions on how to install and use vcf2maf correctly? Thank you very much in advance~~ I am working on macOS Monterey 12.2.1 with an Intel Core i5 processor.

# Answer
> 1 votes
Install it via Anaconda/conda.

```
conda install -c bioconda vcf2maf
conda install -c bioconda/label/cf201901 vcf2maf

```

I personally like Anaconda-navigator for OSX but you would likely go for 'command line version' in part because it will place it in your home directory. Its here https://www.anaconda.com/products/individual and just use the either the above and then its done. The conda environments can be considered once vcf2maf is working.

It is never cool to run Perl on the system command line.

---

Once installed ensure the conda activated

```
conda activate <env name> # if its base I don't think the <> is needed

```

Then simply type

```
vcf2maf.pl

```

**Output**

```
> Usage:
>      perl maf2vcf.pl --help
>      perl maf2vcf.pl --input-maf test.maf --output-dir vcfs

```

Its a Perl script hence it will not recognise maf2vcf without the suffix `pl`. Anaconda doesn't need `perl <perlscript>`.

Oh just to mention I just `conda install vcf2maf` without the channel but I think I've got bioconda on my repo list.

---

I've just uninstalled it. There appears to be an error with the help docs. The help docs (probably PODs) may not have been uploaded and in any case it will run without them. What conda gives you is the full set of dependencies for the script (makes uninstallation a headache) - however you don't need to worry hence I would just use it `maf2vcf.pl --input-maf test.maf --output-dir vcfs` that should work.

---

**Note of explanation: why Anaconda works**

I better explain, the script run by itself is complaining it can't see the dependencies it requires. It doesn't just need them, it needs to see them. Anaconda/conda takes care of all that by downloading the dependencies and ensuring they are all visible. The downside is uninstalling everything afterwards, but that's my problem not the OPs.

---
Tags: vcf, maf
---


title: Multiple Sequence Alignment

text:
# Question
Title: Multiple Sequence Alignment
I applied three multiple sequence alignment methods of muscle, mafft and clistal omega on 30 sequences and 400 sequences . In the case of 30 seq, the output of alignment from the 3 methods was in the same length but in the case of larger data with different length of sequences (input data) the output from alignment methods was different length clustal length is different than the two others. My question is how can i assure that the length from MSA should be equal or different?

# Answer
> 1 votes
The length of the alignment will depend on the model. For example, high gap costs will typically result in shorter alignments, whereas low gap-costs will do the opposite.

Consider the sequences below:

```
AAAAA
TTTTT
GGGGG

```

with a very large gap cost, these sequences will result in the following MSA

```
AAAAA
|||||
TTTTT
|||||
GGGGG

```

but with lower gap costs, we'd observe

```
AAAAA----------
-----TTTTT-----
----------GGGGG

```

Of course, this is an extreme case, but it goes to show that the alignment model will vary how long the resulting alignment is.

# Answer
> 1 votes
It is likely the sequence data in this instance, the more sequences the more chance an issue with the alignment where variations in the gap penalty are observed.

MAFFT outperforms Clustal Omega to some extent over significantly larger alignments \>50 000 sequences from published benchmark studies (FAMSA paper), but this is not relevant here.

---
Tags: variation, multiple-sequence-alignment
---


# Question
Title: tRNAdb insertion symbols. To what they correspond?
I have to manually curate some tRNA sequences in order to have their original sequence without modifications. I am getting these sequences from tRNAdb. In these sequences, I found often "**_**" symbols. They provide this list of symbols and "\_" corresponds to *insertion*. My question is the insertion of what? Are they referring to a reference sequence? In the list they say that in the footnotes they provide more info. However, there are no footnotes. Have you ever had experience with tRNAdb?

# Answer
> 3 votes
I have BLASTed one of the tRNA sequences. Lactobacillus acidophilus, Ile, that has insertion sites.

> GGGCCTATAGCTCAGCTGGTTTAGAGC\_CACGCCTGATAAGCGTG\_GTCGATGGTTCAAGTCCATTTAGGCCCA

It looks like insertion refers to sites that exist in reference DNA but missing from tRNA sequence. But I am not a tRNA expert, it's my interpretation. See the image.

# Answer
> 0 votes
I post the workaround that I eventually found for this problem as it could serve as a reference for the exploration of tRNAdb bioinformatic resource.

If you look at the picture, tRNAdb provides comments about the insertion symbols. Therefore, for the specified position, nucleotides are defined. The original tRNA sequences can be reconstructed from this information. Hope that this problem could be of help to the community!

---
Tags: trna
---


# Question
Title: Taking too long to download metadata file (estimated size 30MB) from Covid-19 Data Portal
I am working on developing a tool to classify covid variants. So I am trying to download SARS-CoV-2 sequence data along with lineage annotation from Covid-19 Data Portal (I could not find a better source that provides both sequence and lineage annotated data.)

Since it is a very large data (given 4,050,433 samples), my idea is to first download the metadata (`.tsv` file) with accession ids and respective lineages from the portal, filter them to keep only the variants I am interested in, and pass the accession\_ids as a text file that I would like to download to CDP File downloader as suggested in the portal.

However, downloading metadata `.tsv` file from the portal is taking forever (download speed is in bytes sometimes). I have tried downloading directly as a file, using `wget` and also using `curl`. I am not sure if there is a way to download it quicker since it is losing connection in between and restarting the download. Just as a side note, the sequence data is downloading pretty fast at 4+Mbps.

I hope my question is clear. Any suggestions/tips on how can I make this happen is very much appreciated. Thanks in advance.

# Answer
> 0 votes
30MB isn't a lot of data at all, given the issue is the EBI servers and not your institutes provision - you would have noticed this before on other projects if it was an institute issue. It therefore suggests someone else is hogging the EBI capacity and there is no provision to bypass this.

The best solution is contact them directly

> ecovid19@ebi.ac.uk

... and state the download size and download speed.

There could be a bot sucking up the database (that is what it sounds like), which would be easy to write given FTP is permitted. They will have the ability to block an IP and NCBI have done this before albeit a long time ago.

Finally filling the user experience form would help EBI which is here. However the issue you described suggests escalation. NCBI are very responsive and I would be surprised if EBI did not respond quickly.

Good luck with the project, sounds interesting.

Generic response below if its just traditional over-demand for the service.

---

I suspect that the beautiful database EBI has provided (and it is amazing) has not anticipated the level of demand for the database size. They are far from being alone.

Modern web architectures, such as Spotify, dedicate enormous effort to matching demand with service provision, notably using Kubernetes. When the server is overloaded they will simply duplicate their core database onto additional machines, which they will be buying in from AWS/GCP. It can even be done dynamically without intervention. NCBI are certainly heading in that direction.

If the EBI wrote this under Kubernetes (it could be a monolith however), the limiting factor is the budget in computer resource allocation (again you can increase Kubernetes provision dynamically). However, if its a monolith this will not be easy and again ultimately someone has pay, i.e. EBI.

Regardless, the fact you can use FTP means EBI are hosting their own service, GCP doesn't permit this and AWS is almost certainly similar. This would mean they have to physically provided additional in-house machines. I visited EBI some time ago and at that time the then director went to some lengths describing how focused EBI were on matching their in-house storage capacity to exponential rise in sequence data. Sudden popularity of a webserver is not trivial issue if it hasn't been anticipated.

---

**Alternatives** Using the NCBI service COVID service is the alternative if there is an urgent need to get the project moving. You may also wish to use GISAID here.

If you want to use the NCBI service, I provided the detailed instructions here: How to programatically download SARS-CoV-2 fasta from NCBI/Genbank via API?

It is however important to understand the possible underlying causes from EBI. I do recall the then EBI director stating their biggest issue was a user simply switching to Genbank, which from their point of view isn't cool given the amount effort going into service provision.

---
Tags: covid-19, data-download, api, curl
---


title: A two color bar plot in R

text:
# Question
Title: A two color bar plot in R
I have raw read counts of miRNAs in two conditions like below

```
> head(long1)
             miRNA Read count Condition
1:  hsa-miR-539-3p         33    living
2: hsa-miR-199b-5p        102    living
3: hsa-miR-146a-3p         27    living
4:   hsa-let-7a-3p        186    living
5:    hsa-miR-1275         54    living
6:  hsa-miR-29a-3p       6464    living
>

```

For each individual miRNAs, I want to show the difference of raw read counts across conditions by

```
p=ggplot(data = long1, aes(x=long1$miRNA, y=long1$`Read count`)) + geom_boxplot(aes(fill=Condition))
p <- p + guides(fill=guide_legend(title="Condition"))
p <- p + xlab("miRNAs") + ylab("Raw mapped reads") + ggtitle("s vs living with cancer")
p + facet_wrap( ~ miRNA, scales="free")+theme_bw() add_pval(p, pairs = list(c("s", "living")), test='wilcox.test')

```

Which returns this

\[!\[enter image description here\]\[1\]\]\[1\]

This plot is not good

How I can a plot like below in which miRNAs are on Y axis and two colours shows the conditions like this

# Answer
> 4 votes
Using Google well is one of the best skillsets you can possibly learn in programming-allied fields. The style of plot you're looking for is a "grouped bar plot". I googled "r grouped bar plot" and found a nice link: https://www.r-graph-gallery.com/48-grouped-barplot-with-ggplot2.html

The key is to use `position=dodge` alongside `color=Condition, x="Read count", y=miRNA` )not sure ow you'd quote the name `"Read count"`, might be better off renaming that column to `readcount`).

---
Tags: r, rna-seq, ggplot2
---


# Question
Title: Is it valid VCF not to 'squash' positions with more than one ALT allele?
I'm seeing output from PBWT that looks like this:

```
3   112841  rs78923776  C   G   .   PASS    RefPanelAF=0.000153988;AN=2;AC=0;INFO=1 GT:ADS:DS:GP    0|0:0,0:0:1,0,0
3   112841  rs78923776  C   T   .   PASS    RefPanelAF=0.078534;AN=2;AC=1;INFO=1    GT:ADS:DS:GP    0|1:0,1:1:0,1,0

```

The first line says that my genome is `CC` at this position, but the second line (for the second alt allele) says that my genome is `CT` at this position. OK, the first line couldn't say this, so it calls me as REF/REF, but this call has to be interpreted in the context of the second call, and can't be taken at face value.

That's why I wonder if this is valid VCF or if it's actually breaking the specification.

The line should be 'squashed' down to the following (I think):

```
3   112841  rs78923776  C   G,T .   PASS    RefPanelAF=0.000153988,0.078534;AN=3;AC=0,1;INFO=1  GT:ADS:DS:GP    0|2:0,0:0:1,0,0

```

# Answer
> 2 votes
It is a valid multi-allelic entry, and most variant callers produce such entries. However, comparison between two such gets more complicated. With biallelic entries, 2 CHR-POS-REF-ALT entries either match or they don't. With multi-allelic entries, they can match completely, partially or not at all, which means you also need to be able to specify how much of a match is a partial match etc. This is a hindrance in annotation processes that annotate using the CHR-POS-REF-ALT fields.

This also introduces cases where an indel can become a multi-allelic entry with a nearby SNV like so:

```
chr1    100    GA    G
chr1    101    A    AA

```

can now be represented as

```
chr1    100    GA    G,AA

```

which can get confusing. Operations are always *more predictable* when VCF has been decomposed and left-aligned.

---
Tags: vcf, file-formats, pbwt
---


# Question
Title: Metagenomics pipeline recommendations for short-read data
*de novo* metagenomics on viral NGS data is a hot-topic. On this site alone at least 4 specific algorithms have been used to identify multi-strain/multi-species for a given data set, however these do not necessarily involve short-read data, and are listed below.

**Questions** What approach is recommended for traditional Illumina multiple infection data both for:

1. for intraspecific (intra-species) data; and
2. for different species co-infecting the patient, but within the same virus family?

Considerations:

* Benchmarking using a common standard (for example)
* Error rates, not miscalling a SNP but producing mosaic artefacts (an artificially mixed genome)
* Suitability of short-read data for intra-species metagenomics.

Programs discussed on Bioinformatics SE:

SAVAGE https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5411778/ Baaijens et al. De novo assembly of viral quasispecies using overlap graphs Genome Research 2017 May; 27: 835–848.

S-aligner (@juanjo75es) https://www.biorxiv.org/content/biorxiv/early/2021/02/02/2021.02.02.429443.full.pdf Juanjo Bermúdez s-aligner: a greedy algorithm for non-greedy de novo genome assembly 2021 bioRχiv February: 1-7

Rey http://denovoassembler.sourceforge.net Sébastien Boisvert et al Genome Biology. 2012 13:R122, December

@LRJoshi here used the Broad Viral NGS pipeline for *de novo* assembly. https://viral-ngs.readthedocs.io/en/latest/assembly.html

---

These are very small genomes \< 30 kb. Quasi-species is a very real phenomenon (loads of mutations), however mixed infections are also well-known. In a mixed infection, producing one genome which is a mix of two genetically distinct viruses is a key concern.

# Answer
Under my experience, the best ones in completeness and N50 are s-aligner, Megahit and rnaSPAdes, but I admit I didn't test the other ones for different reasons. Broad indeed seems to use internally SPAdes and Trinity. I compared the results from these three with reported results in benchmark sets from the other assemblers and the results were better for these three.

From these three I get the best results with s-aligner but it has some issues: it actually requires some practice using it and it is slightly slower for this kind of data (\<30kbp genomes). In the other hand, it usually has lower requirements in memory than other assemblers except maybe Megahit.

Here there are some comparative results (NGA50) using public benchmark sets:

Obviously, I am kinda biased as I am the author of s-aligner but the data doesn't lie.

Getting the best results also may not be the first requirement in some cases, like if you want it fast and with a fast learning curve, and in this case Megahit I think is usually quite good. But in these cases, you also can design a script that will automate s-aligner for your specific needs. You can contact me if you need help.

BTW, s-aligner also gets better results with other kinds of data like bacterial genomes, transcriptomes, or small eukaryotes but always at the cost of requiring more running time or more computational resources (running it in AWS, for example).
> 1 votes

---
Tags: ngs, covid-19, sars-cov-2, metagenome, virus
---


# Question
Title: Bedtools count split reads spanning whole feature
\[Note: this is for bedtools v2.30\]

I want to count the number of RNA-Seq reads that fully cover a given splice junction. For that, I thought of defining a BED feature around the junction, then using `bedtools coverage`.

Here is a simplified situation:

```
1-BASED POS 1|      10|       20|       30|
CHROMOSOME   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
BED FILE              **
BAM READ FWD >>>>>>>>>>
BAM READ REV                  <<<

```

We can generate the corresponding example files:

```
echo -e "I\t9\t11\tmy_position\t.\t+" > single_place.bed

echo -e "@HD\tVN:1.4\tSO:coordinate
@SQ\tSN:I\tLN:30
first_read\t99\tI\t1\t255\t10M\t=\t20\t20\tAATTCCGGAA\tDDDDDDDDDD\tNH:i:1\tHI:i:1\tAS:i:127\tnM:i:0
first_mate\t147\tI\t20\t255\t3M\t=\t1\t20\tATC\tDDD\tNH:i:1\tHI:i:1\tAS:i:127\tnM:i:0" | samtools view -b > fwd_only.bam
samtools index fwd_only.bam

cat single_place.bed
I       9       11      my_position     .       +
samtools view fwd_only.bam
first_read      99      I       1       255     10M     =       20      20      AATTCCGGAA      DDDDDDDDDD      NH:i:1  HI:i:1  AS:i:127        nM:i:0
first_mate      147     I       20      255     3M      =       1       20      ATC     DDD     NH:i:1  HI:i:1  AS:i:127        nM:i:0

```

The per-base count (with the `-d` option) gives the expected result:

```
coverageBed -d -a single_place.bed -b fwd_only.bam
I       9       11      my_position     .       +       1       1
I       9       11      my_position     .       +       2       0

```

And using a filter `-f 1.0` to only keep reads that cover the entire feature indeed removes our read:

```
coverageBed -a single_place.bed -b fwd_only.bam
I       9       11      my_position     .       +       1       1       2       0.5000000
coverageBed -f 1.0 -a single_place.bed -b fwd_only.bam
I       9       11      my_position     .       +       0       0       2       0.0000000

```

However, in practice, I have a lot of spliced reads, so I wish to use the `-split` option. In that case, the filtering doesn't work anymore:

```
coverageBed -split -f 1.0 -a single_place.bed -b fwd_only.bam
I       9       11      my_position     .       +       1       1       2       0.5000000

```

even though the per-base count is still unchanged:

```
coverageBed -d -split -a single_place.bed -b fwd_only.bam
I       9       11      my_position     .       +       1       1
I       9       11      my_position     .       +       2       0

```

Why does `-split` interfere with `-f 1.0`?

### With spliced reads

Note, in practice I'm more interested in cases where the reads can be spliced, such as this:

```
1-BASED POS 1|      10|       20|       30|
CHROMOSOME   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
BED FILE              **
BAM READ FWD >>>>>>>>>>__>>
BAM READ REV                  <<<

```

where `_` is a spliced read. The corresponding BAM can be generated with:

```
echo -e "@HD\tVN:1.4\tSO:coordinate
@SQ\tSN:I\tLN:30
first_read\t99\tI\t1\t255\t10M2N2M\t=\t20\t20\tAAAAAAAAAACC\tDDDDDDDDDDDD\tNH:i:1\tHI:i:1\tAS:i:127\tnM:i:0
first_mate\t147\tI\t20\t255\t3M\t=\t1\t20\tATC\tDDD\tNH:i:1\tHI:i:1\tAS:i:127\tnM:i:0" | samtools view -b > fwd_only.bam
samtools index fwd_only.bam

```

In that situation, it seems I need to use `-split` and `-f 1.0`. I thought of filtering *a posteriori* with e.g. `awk '$10>0.5`, but that doesn't seem to work either: when not using `-d`, the output of `coverageBed` does not separate the reads that cover the feature partially vs fully.

# Answer
> 1 votes
It appears this is an issue that has been raised on Github, and that I had unfortunately missed.

This issue happens in `bedtools coverage` (my examples above), but also in `bedtools intersect`. So, I believe this is actually the same issue as #928. Indeed, switching to bedtools v2.27 solves my issue described above (with v2.30).

Note that this particular issue likely appeared as a consequence of solving issues #773, #750, #673.

My workaround is to use `bedtools intersect` INVERTING the bed and bam, so that I just filter in the reads that match the feature, then count the resulting number of rows.

```
intersectBed -split -bed -wb -F 1.0 \
             -a fwd_only.bam -b single_place.bed | \
  cut -f16 | \
  sort | \
  uniq -c

```

# Answer
> 1 votes
Not quite sure if this is the issue you're having, but I've encountered a similar situation with coverage estimation when trying to map nanopore cDNA reads to the genome. I have a full workflow for mapping oriented long reads here, but here's a summary of the split issue:

> The default options for BEDTools treat sequence deletions (which happen frequently in nanopore reads) as a drop in coverage, which can make exon hunting and coverage calculation more difficult. I submitted a pull request to add an additional ignoreD parameter to the command line to allow cDNA reads with split coverage across introns to ignore deletions when considering coverage; this request has now been incorporated into the main BEDtools repository (as of v2.30.0).

Depending on how your mapping is set up, adding `-ignoreD`to your command line may allow spliced reads to be properly counted for coverage with no deletion gaps:

```
$ bedtools genomecov -bga -strand '+' -split -ignoreD -ibam mapped_BC01.sam
chrX    0       7384244 0
chrX    7384244 7384413 1
chrX    7384413 171031299       0

$ bedtools genomecov -bga -strand '+' -split -ignoreD -ibam mapped_BC01_exon.sam
chr11   0       62551499        0
chr11   62551499        62551595        1
chr11   62551595        62552139        0
chr11   62552139        62553212        1
chr11   62553212        122082543       0

```

---
Tags: bedtools
---


# Question
Title: Potential side effects of replacing read group tags in BAM file
I have a set of BAM files where the read group tags have some (default?) values, i.e.:

```
@RG ID:RG0  LB:LB0  PU:PU0  SM:SM0

```

This creates issues in my downstream analyses, where multiple BAM files with the same SM tag are used.

Samtools provides a command to replace the read group tag. However, I am not sure if there are possible side-effects that I should be aware of, and I might need to remap the BAM file after this replacement. I want to replace `ID` and `SM` with the sample id, which is unique for each BAM file.

Do I need to remap and/or run additional steps, or replacing the RG tags should be sufficient to update the BAM files in a consistent way?

# Answer
> 3 votes
You do not need to remap the files, replacing the read group information with `samtools` is sufficient to deal with this. When you update your pipeline, have it use the sample name to construct the read group information. That way you won't run into this problem again.

# Answer
> 3 votes
I assume you are using GATK software tools for downstream analysis. If that's the case, then you will most likely encounter issues due to missing information on read groups. I believe that GATK can't differentiate between read groups if they are all similar.

In all cases, you can find more information about SAM/BAM and related specifications on the 1.3 The header section as they describe what is mandatory.

So, if you have access to FastQ files, you can add them directly during the alignment process with `bwa mem`:

1. bwa mem

> -R STR Complete read group header line. ’\t’ can be used in STR and will be converted to a TAB in the output SAM. The read group ID will be attached to every read in the output. An example is ’@RG\tID:foo\tSM:bar’. \[null\]

An example command line would be:

```
$bwa mem \
-M \                 // Mark shorter split hits as secondary (for Picard compatibility).
-t "$cpus" \         // number of threads
"$genome" \          // genome reference
"$fastq_r1" \        // R1
"$fastq_r2" \        // provide this line if you have R2
-R "@RG\tID:1\tLB:lib1\tPL:illumina\tSM:$sampleName\tPU:runBarcode" \
> "$sampleName".sam  // Output

```
2. bwa2 mem

> Bwa-mem2 is the next version of the bwa-mem algorithm in bwa. It produces alignment identical to bwa and is ~1.3-3.1x faster depending on the use-case, dataset and the running machine.

The command line is exactly similar to the command line above, except that you would be using `bwa2`.

3. Samtools

As you have mentioned, you can completely overwrite (default mode) the read groups already existing by doing:

```
samtools addreplacerg \
-r "@RG\tID:1\tLB:lib1\tPL:illumina\tSM:$sampleName\tPU:runBarcode" \
-m overwrite_all \         // mode to replace all
-@ 5 \                     // number of cpus
-O b \                     // output bam
-o "$sampleName"_RG.bam \
"$sampleName".bam

```
4. Picard AddOrReplaceReadGroups

```
java -jar $picard AddOrReplaceReadGroups \     // Also available in GATK tools
I="$sampleName".bam \
OUTPUT="$sampleName"_RG.bam \
RGID=1 \
RGLB=library \
RGPL=illumina \
RGPU=runBarcode \
RGSM=${sampleName} \
SORT_ORDER=coordinate

```

---
Tags: sam, samtools
---


# Question
Title: Detect mutation context in a read of a sam file
After sequencing (Illumina) of some DNA, I generated a sam file through alignment of a fastq file (using Bowtie2). Instead of using variant calling programs, I want to know the specific variances, compared to the reference genome, in each specific read.

For example, I'm only interested in C to A mutations. I would like to include reads that contain MORE than 5 C to A mutations. Afterwards I would like to analyze the specific sequence context of that mutation and for doing that I need to know the 3 nucleotides before and after the C to A mutation, which gives a context of NNNCNNN.

I looked at samtools mpileup or bcftools, but I can't find the exact thing that I want to analyze. Any ideas?

Thank you!

# Answer
That kind of custom fiddling with reads and variants is very cumbersome, non-standard and also error-prone. Do a standard variant calling pipeline and then filter for the mutations that you want. Then extract the variant position (so the coordinates) and get the variant context from the reference genome. Using individual reads for this is not advisable as only a properly called variant is reliable. The power of NGS is the validation of an observation by redundancy (=many reads covering a variant). With individual reads your variant can simply be noise, there is hardly a way to tell. Alignment will confidently map the read to the "context" as you call it, so in my opinion there is not really an advantage of calling the "context" directly from individual reads. You could better extract every confident C\>A mutation, then get the three bases up/downstream and then make some kind of position frequency matrix to create a motif.

Example, just for illustration, imagine these were for variants (A in the middle) with three bases up/downstream:

```
CTAATTG
CTCATAG
CTTATAG
CCTAGTG

```

The motif would then be this, you can easily see whether there is some conservation in the "context" of your variant.ad

The variant calling you could do with an end-to-end pipeline such as nf-core sarek, this will eventually give you a VCF file(s) with the variant coordinates. From there simply filter for the variants you want, get the coordinate, add three bases up/downstream to start/end and then use something like `bedtools getfasta` to get the sequences from the reference genome.
> 3 votes

# Answer
If you want to look at a few, specific reads manually, you can do this on VarSome. We have recently released a new feature that allows you to enter single NGS or Sanger reads and get the list of variants detected in that read. For example, here:

In this example, there was one variant identified when comparing to the human hg38 genome assembly:

* CYP24A1(NM\_000782.5):c.1186C\>T which is at chr20:54158136 G⇒A (note that the CYP24A1 gene is on the - strand which is why the C to T variant in HGVS notation is a G to A in the genome). You can then click on the variant link(s) to investigate the variant's effect and annotation.

Obviously, this isn't practical for many reads, but it can be handy if you just want to look at a few manually. More importantly, NGS is very error prone (which is why we sequence at high depth), and you should absolutely not assume that just because a variant is present in a single read it was also actually present in the source DNA you sequenced.

---

Disclaimer: I work for the company that builds VarSome.
> 1 votes

---
Tags: sam, samtools, variant-calling
---


# Question
Title: How to compare differential expression of same genes between two groups?
I need a statistically significant "decision" that DEGs from group A are different (or not) from group B (after drug treatment). For that, I have a mean fold change delta, calculated as the mean of 5 different treated samples minus the mean of the same, but untreated samples.

```
        A.delta    B.delta
gen1    4.9182385  4.5472662
gen2    3.6354761  3.4474844
gen3    3.7592630  3.2838886
gen4    3.9623264  2.8189166
gen5    3.4742171  3.1231050

```

In my vision, there is not much difference in fold change delta between A and B for the `gen1` for example (4.91 is not very far from 4.54). Same for `gen2`. While for `gen4` there might be something interesting. What kind of plot may be the most useful for that kind of comparison?

# Answer
This data is genuinely paired data and you have stated values are a mean. Thus the raw data can be analysed in your previous questions on the same theme described here:

What is the best way to compare treated/control samples?

And here:

What is the best way to compare differences in DEG expression between two populations?

If this is 'ANOVAerable' (its not a word, but if ANOVA is appropriate) than a two-way ANOVA would provide all the answers (however see below regarding the trend in the data). This analysis will perform a 'A' versus 'B' and each 'gen' compared to the other, as well as the interaction between A/B and 'gen'. You would need the raw data, you can't do anything with means (unless you are assuming the distribution). The distribution within each 'gen' is important (if its normal - conforms to the normal distribution, you're in business).

There is nothing further to add and I do feel we need to draw this series of questions to a close.

A regression model could be explored because there is a clear trend in both data sets, but I'm uncertain it will be significant. Within ANOVA however this might trap the outlier, but that is more complicated (you standardise against the trend thus 'gen1' is comparable to 'gen5' and therein 'gen4' can be identified as an outlier). This was touched on in the answer to the previous question (ols and python)
> 0 votes

---
Tags: differential-expression
---


# Question
Title: What does it mean if aligned genomes are not of the same length?
I am aligning 4 genomes of *Psuedomonas chlororaphis* using progressive mauve software. Command used:

```
progressiveMauve --output=outputfile --disable-backbone chl189.fas chl17415.fas chl50083.fas chlTAMOak81.fas

```

At the end, the resulting aligned genomes are not of the same length. What could be the potential problems?

This is a global alignment. Which is why im surprised that they aren't all the same length after alignment? Also, while analysing the subsequent genomes, I noticed that some region which were present in only 2/4 genomes were misaligned in the end result. How do I proceed troubleshooting this problem?

I am a student and new to this field, could you please tell me how to proceed with manually checking to see whats wrong?

# Answer
> 1 votes
I’ve seen cases like this. For me, there were plasmids at end that didn’t get aligned.

---
Tags: sequence-alignment, genome, mauve
---


# Question
Title: Database to download with information of bacteria?
I am looking for tidy tabular dataset that contains known bacteria and fungi down to the species level. Basically a table that includes genus and species names and other higher-level information. For example, whether it is gram-positive or negative, or anaerobic etc. A database with API would work too.

I found https://bacteria.ensembl.org/, but this is not so easy to use/automate.

# Answer
You can try Entrez in Biopython. Although, there is hassle with additional information.

```
from Bio import Entrez

Entrez.email = "your.letterbox@any.com"
Entrez.api_key = "your_key"

TaxID = Entrez.read(Entrez.esearch(db="taxonomy", term="Hepatovirus A"))["IdList"][0]
Entrez.read(Entrez.efetch(db="taxonomy", id=TaxID))

```
> 1 votes

---
Tags: python, database, public-dataset
---


title: Installing modules in perl for Circos

text:
# Question
Title: Installing modules in perl for Circos
Im trying to install the missing modules of Circos without success.

The error code is:

```
cpan[1]> install Config::General
Reading '/home/asg28/.cpan/Metadata'
  Database was generated on Mon, 14 Mar 2022 03:29:03 GMT
Running install for module 'Config::General'
Fetching with HTTP::Tiny:
http://www.cpan.org/authors/id/T/TL/TLINDEN/Config-General-2.63.tar.gz
Fetching with HTTP::Tiny:
http://www.cpan.org/authors/id/T/TL/TLINDEN/CHECKSUMS
Checksum for /home/asg28/.cpan/sources/authors/id/T/TL/TLINDEN/Config-General-2.63.tar.gz ok
Scanning cache /home/asg28/.cpan/build for sizes
............................................................................DONE
'YAML' not installed, will not store persistent state
Configuring T/TL/TLINDEN/Config-General-2.63.tar.gz with Makefile.PL
Checking if your kit is complete...
Looks good
Generating a Unix-style Makefile
Writing Makefile for Config::General
Writing MYMETA.yml and MYMETA.json
  TLINDEN/Config-General-2.63.tar.gz
  /usr/bin/perl Makefile.PL INSTALLDIRS=site -- OK
Running make for T/TL/TLINDEN/Config-General-2.63.tar.gz
Can't exec "make": No such file or directory at /usr/share/perl/5.30/CPAN/Distribution.pm line 2225.
  TLINDEN/Config-General-2.63.tar.gz
  make -- NOT OK
  No such file or directory
Failed during this command:
 TLINDEN/Config-General-2.63.tar.gz           : make NO

```

I'm running Ubuntu in Windows and already tried in another computer

# Answer
> 1 votes
The simplest solution is to install conda/Anaconda

Then type any one of the following. Start with the top command and if it is unavailable go to the next commend down - only needs work once:

```
conda install -c bioconda circos
conda install -c bioconda/label/broken circos
conda install -c bioconda/label/cf201901 circos

```

This is a preferred solution because the `make` will already be configured along with any other dependency and in addition you can run it directly on Windows without an intervening simulator as described here: https://docs.anaconda.com/anaconda/install/windows/

You could install `brew` but this would be a bit messy,

```
/bin/bash -c "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/HEAD/install.sh)"

```

because you would be using a perl environment within a simulator.

The alternative is to do what @user324810 and install `make` explicitly and is described here, https://stackoverflow.com/questions/23293589/perl-module-install-error-cpan

The error is easy enough but there will be future dependencies which will not be installed and a once an for all solution is a better option and diagnosing and installing each dependency in turn.

---
Tags: perl, molecular-genetics, circos
---


# Question
Title: On the same strand, for the same gene, can exons be overlapping?
I want to get a set of exon regions for each protein coding gene.

I extracted a set of relevant information (chromosome, start, end, gene ID, gene name, gene type, exon number and exon ID) from a GTF using:

```
awk '$3 == "exon"' annotation.gtf \
  | gtf2bed \
  | cut -f1-4,10 \
  | awk 'match($0, /gene_name \"[a-zA-Z0-9\_\-.\/]*\"/, gene_name)
    && match($0, /gene_type \"[^"]*\"/, gene_type)
    && match($0, /exon_number [0-9]*/, exon_number)
    && match($0, /exon_id \"[^"]*\"/, exon_id) {
      print $1,$2,$3,$4,gene_name[0],gene_type[0],exon_number[0],exon_id[0];
    }' \
  | tr -d '"' \
  | tr ' ' $'\t' \
  | cut -f1-4,6,8,10,12 \
  | grep 'protein_coding' \
  > annotation.exons.bed

```

However, during a sanity check, I noticed that for some genes there are some exons with overlapping coordinates, same gene information and exon number, but different exon ID, e.g.:

```
chr2    9423650 9424052 ENSG00000119203.14      CPSF3   protein_coding  1       ENSE00001829522.1
chr2    9423653 9423823 ENSG00000119203.14      CPSF3   protein_coding  1       ENSE00001830696.2
chr2    9423883 9424960 ENSG00000119203.14      CPSF3   protein_coding  1       ENSE00001932298.1

```

I did not directly work with exons in the past, but I am aware of alternative splicing. My understanding was that exons can be skipped during alternative splicing, but I assumed that exon number was unique for bases, within the same gene (on the same strand).

In the example above, the first 2 listed exons have the same exon number, overlapping (but not exactly the same) bases, and different exon ID. Is this due to some additional information that I can retrieve from the GTF, i.e., from:

```
chr2    9423650 9424052 ENSG00000119203.14      .       +       HAVANA  exon    .       gene_id "ENSG00000119203.14"; transcript_id "ENST00000475482.5"; gene_type "protein_coding"; gene_name "CPSF3"; transcript_type "protein_coding"; transcript_name "CPSF3-203"; exon_number 1; exon_id "ENSE00001829522.1"; level 2; protein_id "ENSP00000419744.1"; transcript_support_level "4"; hgnc_id "HGNC:2326"; tag "mRNA_end_NF"; tag "cds_end_NF"; havana_gene "OTTHUMG00000090415.5"; havana_transcript "OTTHUMT00000353331.3";
chr2    9423653 9423823 ENSG00000119203.14      .       +       HAVANA  exon    .       gene_id "ENSG00000119203.14"; transcript_id "ENST00000238112.8"; gene_type "protein_coding"; gene_name "CPSF3"; transcript_type "protein_coding"; transcript_name "CPSF3-201"; exon_number 1; exon_id "ENSE00001830696.2"; level 2; protein_id "ENSP00000238112.3"; transcript_support_level "1"; hgnc_id "HGNC:2326"; tag "basic"; tag "MANE_Select"; tag "appris_principal_1"; tag "CCDS"; ccdsid "CCDS1664.1"; havana_gene "OTTHUMG00000090415.5"; havana_transcript "OTTHUMT00000206843.2";

```

?

Is this because the different exons belong to different transcripts? However, shouldn't different exons still be non-overlapping one against each other?

# Answer
> 5 votes
I'm not quite sure what you're asking, but pretty much anything can happen in terms of exons and transcripts:

* Two different genes encoded on different strands can share the same genomic region as an exon
* Two different genes encoded on *the same strand* can share the same genomic region as an exon (including identical exon boundaries)
* A gene can have multiple transcript isoforms that share the same genomic region on the same frame
* A gene can have multiple isoforms that share the same genomic region in different frames
* Genes can have variable regions which are modified during somatic development (e.g. T-cell and B-cell receptors)
* Gene transcripts can be spliced together in different orders in different isoforms
* Exons from one transcript isoform can overlap with exons from another isoform (but not necessarily cover identical regions)
* Stop codons can occasionally be "read-through", resulting in a longer protein/peptide than the expected translation

As an example of different isoforms within the same gene exhibiting different translation mechanisms, here are graphical representations of annotated transcript isoforms for Gapdh on the mouse GRCm38 assembly, created in JBrowse. I wouldn't put complete trust in any of these (especially not the ones that have no defined intronic regions), but it's a useful demonstration of the types of things that can happen. Compare, Gapdh-207, where an untranslated region starts halfway through an exon, and Gapdh-209, where the translated protein starts halfway through an untranslated exon in Gapdh-207:

If there is a situation where it seems like exons *from the same isoform* have overlapping regions, it might be worth looking elsewhere in the genome to see if there's a transcribed gene copy somewhere. That, at least, doesn't seem to be something I'd expect to happen through the normal translation process.

Also, don't expect annotations to be complete. Even on something as \[presumably\] well-studied as the mouse mitochondrial genome, we can identify clear transcript start/stop points from long polyA cDNA sequences that are completely absent in the genome annotations. This is most obvious for the completely unannotated reverse strand (e.g. transcript starting at around 16,200bp), but also happens on the forward strand as well (e.g. coverage steps within mt-Nd5):

---
Tags: proteins, sequence-annotation, gtf, rna-splicing, exon
---


# Question
Title: Nextflow: cannot get file pairs from Channel.fromFilePairs
I'm remarkably still struggling with some basic `Nextflow` concepts. I'm simply trying to collect all pairs of `.cram` and `.crai` files from a directory using `Channel.fromFilePairs(...)`, but I get no output when I try to `view()` the `Channel`.

**Question:** Why does my `Channel.fromFilePairs(...)` command not return file pairs, but I get a list of files when using `Channel.fromPath(...)`?

As far as I can tell, I'm following the same pattern from the following examples:

1. Nextflow fromFilePairs
2. .bam/.bai example (except without the extra code)

### Code

```
nextflow.enable.dsl=2

workflow {

    /*
     * Channel.fromPath returns the expected files and prints them to
     * terminal.
     */
    Channel.fromPath('/path/to/adsp/cram/snd10000/A-CUHS-CU000*.cram')
        | view()

    /*
     * Channel.fromFilePairs does not print anything to terminal.
     */
    Channel.fromFilePairs('/path/to/adsp/cram/snd10000/A-CUHS-CU000*.{cram,crai}')
        | view()
}

```

### Output (note, no file tuples/pairs)

```
N E X T F L O W  ~  version 21.04.3
Launching `pairs.nf` [pensive_plateau] - revision: f917a5c41b
/path/to/adsp/cram/snd10000/A-CUHS-CU000563-BL-COL-47451BL1_vcpa1.1.cram.crai
/path/to/adsp/cram/snd10000/A-CUHS-CU000387-BL-COL-38200BL2_vcpa1.1.cram.crai
/path/to/adsp/cram/snd10000/A-CUHS-CU000018-BL-COL-46250BL1_vcpa1.1.cram.crai
/path/to/adsp/cram/snd10000/A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1.cram
/path/to/adsp/cram/snd10000/A-CUHS-CU000563-BL-COL-47451BL1_vcpa1.1.cram
/path/to/adsp/cram/snd10000/A-CUHS-CU000387-BL-COL-38200BL2_vcpa1.1.cram
/path/to/adsp/cram/snd10000/A-CUHS-CU000194-BL-COL-38470BL1_vcpa1.1.cram
/path/to/adsp/cram/snd10000/A-CUHS-CU000194-BL-COL-38470BL1_vcpa1.1.cram.crai
/path/to/adsp/cram/snd10000/A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1.cram.crai
/path/to/adsp/cram/snd10000/A-CUHS-CU000018-BL-COL-46250BL1_vcpa1.1.cram

```

# Answer
> 3 votes
**TLDR**, use:

```
nextflow.enable.dsl=2

workflow {

    Channel.fromFilePairs('/path/to/adsp/cram/snd10000/A-CUHS-CU000*.cram{,.crai}')
        | view()
}

```

Results:

```
N E X T F L O W  ~  version 20.10.0
Launching `xscript.nf` [special_caravaggio] - revision: a4e53657a0
[A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1, [/path/to/adsp/cram/snd10000/A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1.cram, /path/to/adsp/cram/snd10000/A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1.cram.crai]]
[A-CUHS-CU000194-BL-COL-38470BL1_vcpa1.1, [/path/to/adsp/cram/snd10000/A-CUHS-CU000194-BL-COL-38470BL1_vcpa1.1.cram, /path/to/adsp/cram/snd10000/A-CUHS-CU000194-BL-COL-38470BL1_vcpa1.1.cram.crai]]
[A-CUHS-CU000018-BL-COL-46250BL1_vcpa1.1, [/path/to/adsp/cram/snd10000/A-CUHS-CU000018-BL-COL-46250BL1_vcpa1.1.cram, /path/to/adsp/cram/snd10000/A-CUHS-CU000018-BL-COL-46250BL1_vcpa1.1.cram.crai]]
[A-CUHS-CU000387-BL-COL-38200BL2_vcpa1.1, [/path/to/adsp/cram/snd10000/A-CUHS-CU000387-BL-COL-38200BL2_vcpa1.1.cram, /path/to/adsp/cram/snd10000/A-CUHS-CU000387-BL-COL-38200BL2_vcpa1.1.cram.crai]]
[A-CUHS-CU000563-BL-COL-47451BL1_vcpa1.1, [/path/to/adsp/cram/snd10000/A-CUHS-CU000563-BL-COL-47451BL1_vcpa1.1.cram, /path/to/adsp/cram/snd10000/A-CUHS-CU000563-BL-COL-47451BL1_vcpa1.1.cram.crai]]

```

---

---

When using the fromFilePairs factory method, the trick is to ensure that the pattern you provide can return a *common* prefix in order for the files to be grouped as expected. An easy way to check to see if this is working as expected is to directly call the readPrefix helper function. For example:

```
test_cram = file('/path/to/adsp/cram/snd10000/A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1.cram')
test_crai = file('/path/to/adsp/cram/snd10000/A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1.cram.crai')

pattern = '/path/to/adsp/cram/snd10000/A-CUHS-CU000*.{cram,crai}'

println("My CRAM prefix is: ${Channel.readPrefix(test_cram, pattern)}")
println("My CRAI prefix is: ${Channel.readPrefix(test_crai, pattern)}")

```

You can see that these results do *not* produce a common prefix:

```
N E X T F L O W  ~  version 20.10.0
Launching `script.nf` [lonely_jang] - revision: 428df7e55e
My CRAM prefix is: A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1
My CRAI prefix is: A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1.cram

```

But if you change your pattern definition to:

```
pattern = '/path/to/adsp/cram/snd10000/A-CUHS-CU000*.cram{,.crai}'

```

We can check that it returns a common prefix:

```
N E X T F L O W  ~  version 20.10.0
Launching `script.nf` [maniac_volta] - revision: a8acf7b03f
My CRAM prefix is: A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1
My CRAI prefix is: A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1

```

---
Tags: nextflow
---


title: How to use GEO_SPHERE?

text:
# Question
Title: How to use GEO_SPHERE?
I have some problems with geosphere package which extends/applies the Beast Bayesian MCMC package for calculating phylogenetic trees. Please, could anyone give any hints about at least one of them?

Firstly, I cannot reproduce the tutorial. I followed the instructions step-by-step but my results are too different (Fig 1) from specified in it. What could I do wrong? Secondly, it is oft occurred that my mean and standard deviation increase too high (Fig 2). I am aware that this is bad but what is the right way not to allow such situations? Just put an upper border? In this case, what is the reasonable border? Thirdly, why is a site model of the geodata partition set as the first option in a list of models for amino acids or nucleotides substitution, & I cannot pick any other (Fig 3 & 4)? What do precision and fast checkbox mean? Is there any way to get under the hood of model to vary more parameters? Finally, are there any documentation, guides, tutorials &c. about geosphere besides the aforementioned tutorial? I did not find any on its GitHub repository.

Fig 1

Fig 2

Fig 3

Fig 4

My trees:

# Answer
As it is turned out from the author's words:

1. This is a normal result for the tutorial.
2. It is recommended to border the deviation of geomodel to 3 powers.
3. There is set up a homogenous diffusion model but BEAUTi does not show this correctly. Precision is defy τ as bt^-1 where t – a branch length, b – precision. It can be an arbitrary number that the quotient of it divided by the potential tree height would much lower than 1. Otherwise, the next approximation has some deterioration:

Threshold is an angle (in degrees 0...180) below which there is no contribution to the probability density. The fast check box turns up approximation of arccos of angles close to zero (cos(x) \> 0.9). The error is at the most 1e-10. More hidden parameters can be found through XML-editors in prepared files.

By the way, the author, Remco Bouckaert, has released new version recently where we can specify regions for clades explicitly in BEAUTi.

I hope that this information will be useful for anyone like me who is stuck in GEO\_SPHERE.
> 1 votes

# Answer
This may not be the exact answer but its a pretty good guess and suggests an absence of `getters` and `setters` for this particular issue. Thus I accept the error does not look informative.

I know 100% what is wrong with the MCMC tree model you performed on HBV, hepatitis B virus.

The `ESS` are mostly \<\<200, that is where they are in red. The MCMC needs to run for longer or a simpler model (probably both). For example, the `prior` is 4 which is way too low.

In particular `sprecision.geo` is `-` and that does not look good for the calculation you are performing because it suggests it (the `geo` object is not present and looks remarkably like the error you are receiving). Simply this will be written in Java and that value it will see as an object. Whether this is singly due to insufficient runtime, I don't know, but the ESS values are a pretty good guess at diagnosing the situation.

What might be happening is that if the values in red are too low and insufficient to calculate `sprecision.geo` hence its `-` and the absence of this value is preventing the `diffusionModel.s:geo` from being calculated.

Do consider Beast is a hugely impressive but complicated population dynamics Bayesian simulation which is further complicated because you are stacking a further level of parameters onto it.

---

I have responded to the comments below here:

HBV is known to be a difficult data set in phylogeny, it has a very weird genome and if further iterations are resulting in errors with the basic MCMC something is definitely not right with the core calculation. 'Too much MCMC' isn't a valid error in a Bayesian calculation. I would simply switch data set and try again. The ESS is simply too low for any further analysis and that's why its in red.

One of the photos provided doesn't seem to stack up.
> 0 votes

---
Tags: phylogenetics, phylogeny, virus, beast
---


title: Can't load files into PopART

text:
# Question
Title: Can't load files into PopART
I am trying to create haplotype networks in PopART (Population genetics with Reticulated Trees\*) which I have successfully done many times in the past.

Two weeks ago the program was working fine for me but recently I am unable to load files. What happens is that when I go to load a nexus file through the "nex" button, or through the "file/import/alignment" or "file/open" options nothing loads.

There is no error message, instead nothing happens. This occurs with both files I have recently created and with files that I have successfully loaded in the past. I've sent files to other people and everyone has been able to load the files using PopART on their computers. I am running macOS 10.15.6 and have tried turning my computer off and on as well as uninstalling and reinstalling PopART. I also thought that maybe Java wasn't up to date on my computer but I don't think PopArt runs with Java. Any help would be very appreciated!

---

PopArt builds reticulated HAplotype networks such shown below

# Answer
> 0 votes
*Answer from @user438383, converted from comment:*

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

# Answer
> 0 votes
PopArt have a mailing list popartusers+subscribe@googlegroups.com

If you go here you can subscribe to it. This sounds like an issue with Mac compatibility. OS 10.15.6 isn't a version I understand and suspect there's a typo. Currently I believe it is OSX 12.1, thus OS10.12.1 PopArt was written for OSX10.6 and upwards. HOWEVER, there can be issues with major system upgrades so Monterey OSX12, I think was a major revision. PopArt is no longer supported so it would be sensitive to changes in Mac OS.

My suggestion is to try this on a Windows machine, if you get the same error the suspicion is there has been a small change to the input format.

---
Tags: phylogenetics, phylogeny, haplotypenetwork
---


# Question
Title: Microbial diversity analysis using whole-genome metagenomic data
I have data, obtained from a single metagenomic DNA sample, that consists of two MiSeq FASTQ files (R1 and R2) that I merged using PEAR.

Now I want to estimate the abundances of the bacteria taxa to generate a figure like this one:

Figure from: *Panosyan, Hovik, and Nils‐Kåre Birkeland. "Microbial diversity in an Armenian geothermal spring assessed by molecular and culture‐based methods." Journal of basic microbiology 54.11 (2014): 1240-1250.*

The problem is that there wasn't a step of amplification of the 16S region as the goal of the sequencing was to discover new genes. I've already isolated 16S reads from my sample using SortMeRNA, but it seems like softwares that do OTU picking, taxonomic assignment and diversity analyses (such as mothur and QIIME) require that all the reads come from the same region of the 16S gene.

Is there a way of using these 16S reads that I've filtered using SortMeRNA in a diversity analysis using mothur/QIIME?

# Answer
New methods have been published in the time since this question was asked. Now, the relative abundance of organisms across different taxonomic levels can be estimated directly from metagenomic shotgun sequencing data without the need to identify and pull rRNA genes. Here are a few examples of popular tools:

* Similar to the approach described in the question, **mOTUs** employs a database of single-copy phylogenetic marker genes to create taxonomic profiles from short-read metagenomic data.
* **Kraken** associates short genomic substrings (*k*-mers) with the lowest common ancestor in its database. The output can then be processed with **Bracken** to get estimates of relative abundance at any taxonomic level, and the Bracken output table can then be passed to a tool like **bracken\_plot** to visualize differences in relative abundance between samples.
* Like mOTUs, **MetaPhlAn** uses a large marker gene database to assign reads to bacterial, archaeal, viral, and eukaryotic clades. This allows for relative abundance estimation in terms of number of cells. If their paper is to be believed, MetaPhlAn3 outperforms both mOTUs2 and Kraken2/Bracken in terms of abundance profile accuracy.

> In addition to more accurate species detection, MetaPhlAn3 also quantified taxonomic abundance profiles more accurately compared to MetaPhlAn 2, mOTUs2, and Bracken based on Bray-Curtis similarity in most datasets ... While it was slightly outperformed by mOTUs2 (only in high-recall mode) on the synthetic mouse gut dataset, even in this case, correlation-based measures (Pearson Correlation Coefficient between estimated and expected relative abundances) found MetaPhlAn3 to be more accurate (r = 0.73) than the other considered profilers (MetaPhlAn2 r = 0.63, mOTUs2 precision r = 0.60, mOTUs2 recall r = 0.71, Bracken r = 0.43).

Although it is now dated, *How to Choose Your Metagenomics Classification Tool* from the Salzberg lab lists more options and links the associated publications, including **Centrifuge** and **Pavian**, the latter of which is geared specifically toward visualization of metagenomic classification data.
> 1 votes

# Answer
This is maybe a very simplified approach but have you thought about simply blasting them against the rRNA/ITS database of NCBI/ SILVA with ie. diamond?

First download the database files from either of the database:

SILVA: https://www.arb-silva.de/no\_cache/download/archive/current/Exports/ NCBI: https://ftp.ncbi.nlm.nih.gov/blast/db/

Second: Use DIAMOND to align fastas of interest (in your case these will be the 16S rRNA Sequences. These will give you hits per fasta which you can save into tsv or csv format and parse in excel.

To get the abundances, take your fastas of interest and align them back to your raw files (paired reads). You can use BOWTIE + SAMTOOLs for this.

Concatenate the results from the DIAMOND output with the abundances and voila, you should get a nice table with taxa and their abundances. If you have more samples and you would like to compare them, ie. compare the taxa composition. I would strongly suggest that you do some further normalisation ie. normalise the abundances with sequencing depths.

Maybe it's a long way around it but in theory, it should do the job.
> 0 votes

---
Tags: metagenome, taxonomy
---


# Question
Title: LD Clumping vs. Conditional Analysis for independent associations
I am working with some GWAS summary statistics. My ultimate goal is to extract variants to use as instrument variables for Mendelian Randomization analysis, which requires me to find SNPs that are independently associated with the exposure variable of interest.

My current workflow involves extracting SNPs from summary statistics with P-value \< 5E-08, then performing LD Clumping with a stringent R^2 cutoff (R^2\<0.05) within a locus. Is this sufficient to ensure that my variants are independently associated, or should I also plan on running conditional analysis(e.g. GCTA-COJO) in addition to/in lieu of clumping for my variants?

# Answer
LD Clumping with $r^2$ filtering seems like a good start, but bear in mind that population structure has a significant effect on association statistics.

I recommend carrying out separate PCAs (or some other dimensional reduction) on your variable groups (e.g. risk/non-risk, or high/medium/low) to make sure that their genomic population structure is similar. Alternatively / equivalently, carry out a whole-group reduction, and confirm that the coordinate centres for each variable group are similar.

If the population structures are *not* similar (which is very likely), then it's important to \[also\] exclude structure-specific SNPs from the analysis, as these will confound the results of your analysis. Use a clustering method to identify people with similar genetic structure, and treat the different clusters as different trait groupings for identifying those SNPs.
> 1 votes

---
Tags: gwas
---


# Question
Title: Convert FASTA to FASTQ with dummy quality scores
I have a FASTA file which I would like to convert into FASTQ format as the tool I want to use my data in requires it in FASTQ format. So dummy quality scores are fine.

Note: I am not using any functionality of the tool that will require the quality score.

# Answer
> 4 votes
One way of doing this is with two subcommands from the `pyfastaq` suite.

```
fastaq to_fake_qual in.fasta - | fastaq fasta_to_fastq in.fasta - out.fastq

```

The first tool, `to_fake_qual`, creates fake quality scores (default 40) for each base and the `-` sends that file (`.qual`) to stdout. The second tool, `fasta_to_fastq`, consumes both the original fasta and the quality scores coming from stdin and turns these into a fastq file.

# Answer
> 6 votes
A simple Biopython solution:

```
from Bio import SeqIO

for r in SeqIO.parse("myfile.fa", "fasta"):
    r.letter_annotations["solexa_quality"] = [40] * len(r)
    print(r.format("fastq"), end='')

```

Example:

```
$ cat << EOF > myfile.fa
>1
ATG
>2
ATGTAGA
EOF
$ python3 myscript.py
@1
ATG
+
III
@2
ATGTAGA
+
IIIIIII

```

# Answer
> 2 votes
I think the python package bioconvert does what you want if you don't provide qualities: https://github.com/bioconvert/bioconvert/blob/master/bioconvert/fasta2fastq.py

The following outputs lots of warnings for me, but does generate a fastq file `test.fq` from an existing fasta file `test.fa` (using dummy quality "I", and adding "None" to sequence headers that do not have a comment part in the original fasta):

```
bioconvert fasta2fastq test.fa test.fq

```

(This requires installing bioconvert using `pip` beforehand: `pip install bioconvert`.)

# Answer
> 2 votes
Here's a quick one-liner with no dependencies other than the UNIX core utils and Perl.

```
cat seqs.fasta | paste - - | perl -ne 'chomp; s/^>/@/; @v = split /\t/; printf("%s\n%s\n+\n%s\n", $v[0], $v[1], "B"x length($v[1]))' > seqs.fastq

```

Caveat: assumes your Fasta sequences are not wrapped, i.e. two lines per record.

# Answer
> 0 votes
I have just uploaded a Python tool to convert fasta to fastq: https://github.com/Mostafa-MR/A\_tool\_to\_convert\_fasta\_to\_fastq

Please note this is presented directly from a Jupyter notebook, i.e. it is a ipynb.

---
Tags: fasta, fastq
---


# Question
Title: Is it possible to reconstruct an MSA with PSSM?
I am thinking of reconstructing an multiple sequence alignment (MSA) of protein sequences from a position-Specific Scoring Matrices (PSSM).

Is it possible? I suppose co-evolution information is lost in PSSM so maybe it is difficult?

# Answer
In short I'm afraid it isn't possible because you lose the linkage per Haplotype which is the core information and vital for any population genetics or phylogenetic analysis.

Position-specific scoring matrix e.g. A .60 C .20 G .20 T 0 I tried to produce a position scoring matrix, but it didn't work out, anyway the allele frequency for each e.g. amino acid are given for each alignment position.

You could produce a consensus sequence which, albeit be hypothetical, it might be of some use to understand the generic behaviour of all the sequences in your sample.

Here's a cartoon describing linkage and its importance,

For a multiple sequence alignment of course AB and ab are amino acid residues within the same locus. PSSM could be used to construct a matrix to analyses an established alignment, but that's a fairly complicated Phylogenetics calculation (and a separate question).
> 1 votes

---
Tags: phylogenetics, proteins, phylogeny, pssm
---


# Question
Title: Installing snakemake following official docs just won't work. Help needed
Following the docs just won't work for me.

```
conda install -n base -c conda-forge mamba
conda activate base
mamba create -c conda-forge -c bioconda -n snakemake snakemake

```

This is the Dockerfile I use:

```
FROM ubuntu:20.04
ENV PATH /opt/conda/bin:${PATH}
ENV LANG C.UTF-8
ENV SHELL /bin/bash
COPY aws/docker/Miniconda3-py38_4.10.3-Linux-x86_64.sh /tmp/
RUN bash /tmp/Miniconda3-py38_4.10.3-Linux-x86_64.sh -b -p /opt/conda
RUN conda install -y -c conda-forge mamba
RUN mamba create -q -y -c conda-forge -c bioconda -n snakemake snakemake snakemake-minimal --only-deps

```

The error is:

```
WARNING Extracted package cache '/opt/conda/pkgs/_openmp_mutex-4.5-1_gnu' has invalid size
WARNING Extracted package cache '/opt/conda/pkgs/_openmp_mutex-4.5-1_gnu' has invalid SHA-256 checksum
WARNING Package tarball '/opt/conda/pkgs/_openmp_mutex-4.5-1_gnu.tar.bz2' is invalid
  Package                              Version  Build                Channel                     Size
───────────────────────────────────────────────────────────────────────────────────────────────────────
  Install:
───────────────────────────────────────────────────────────────────────────────────────────────────────

  ...Package list omitted

  Ignored:
───────────────────────────────────────────────────────────────────────────────────────────────────────

  = snakemake                            6.9.1  hdfd78af_0           bioconda/noarch          Ignored
  = snakemake-minimal                    6.9.1  pyhdfd78af_0         bioconda/noarch          Ignored

  Summary:

  Install: 215 packages
  Ignored: 2 packages

  Total download: 256 MB

───────────────────────────────────────────────────────────────────────────────────────────────────────

ERROR   File not valid: SHA256 sum doesn't match expectation "/opt/conda/pkgs/gettext-0.19.8.1-h73d1719_1008.tar.bz2"
        Expected: 1bb53c99b4943d210c881aad9158fb0235b348498bad1a7076d1f2bef6671922
        Actual: e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855

ERROR   File not valid: SHA256 sum doesn't match expectation "/opt/conda/pkgs/gtk2-2.24.33-h539f30e_1.tar.bz2"
        Expected: 2a0aef09a947e3768a426b0ea414ddd1e015c9804a46de619a5e67a689ea1dec
        Actual: e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855

ERROR   File not valid: SHA256 sum doesn't match expectation "/opt/conda/pkgs/botocore-1.21.64-pyhd8ed1ab_0.tar.bz2"
        Expected: ee5501c74eefd0396f141593933ef2cc3419ad21624ed140bbc1e04944af6e31
        Actual: e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855

ERROR   File not valid: SHA256 sum doesn't match expectation "/opt/conda/pkgs/libopenblas-0.3.18-pthreads_h8fe5266_0.tar.bz2"
        Expected: 6d67f0f7bda1503af83659ca45d8c8ba34d94091fa4cd5b1ecd85332c14db3cc
        Actual: e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855

ERROR   File not valid: SHA256 sum doesn't match expectation "/opt/conda/pkgs/scipy-1.7.1-py39hee8e79c_0.tar.bz2"
        Expected: 536500b79a5ccc91f6059fcaa6d7e959bddb37a0ed7a930797375c1ed51a6af1
        Actual: e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855

# >>>>>>>>>>>>>>>>>>>>>> ERROR REPORT <<<<<<<<<<<<<<<<<<<<<<

    Traceback (most recent call last):
      File "/opt/conda/lib/python3.8/site-packages/conda/exceptions.py", line 1079, in __call__
        return func(*args, **kwargs)
      File "/opt/conda/lib/python3.8/site-packages/mamba/mamba.py", line 926, in exception_converter
        raise e
      File "/opt/conda/lib/python3.8/site-packages/mamba/mamba.py", line 920, in exception_converter
        exit_code = _wrapped_main(*args, **kwargs)
      File "/opt/conda/lib/python3.8/site-packages/mamba/mamba.py", line 878, in _wrapped_main
        result = do_call(args, p)
      File "/opt/conda/lib/python3.8/site-packages/mamba/mamba.py", line 757, in do_call
        exit_code = create(args, parser)
      File "/opt/conda/lib/python3.8/site-packages/mamba/mamba.py", line 660, in create
        return install(args, parser, "create")
      File "/opt/conda/lib/python3.8/site-packages/mamba/mamba.py", line 609, in install
        downloaded = transaction.prompt(repos)
    RuntimeError: Found incorrect download: libopenblas. Aborting

`$ /opt/conda/bin/mamba create -q -y -c conda-forge -c bioconda -n snakemake snakemake snakemake-minimal --only-deps`

  environment variables:
                 CIO_TEST=<not set>
               CONDA_ROOT=/opt/conda
           CURL_CA_BUNDLE=<not set>
                     PATH=/opt/conda/bin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin
                          :/bin
       REQUESTS_CA_BUNDLE=<not set>
            SSL_CERT_FILE=<not set>

     active environment : None
       user config file : /root/.condarc
 populated config files :
          conda version : 4.10.3
    conda-build version : not installed
         python version : 3.8.10.final.0
       virtual packages : __linux=5.10.16.3=0
                          __glibc=2.31=0
                          __unix=0=0
                          __archspec=1=x86_64
       base environment : /opt/conda  (writable)
      conda av data dir : /opt/conda/etc/conda
  conda av metadata url : None
           channel URLs : https://conda.anaconda.org/conda-forge/linux-64
                          https://conda.anaconda.org/conda-forge/noarch
                          https://conda.anaconda.org/bioconda/linux-64
                          https://conda.anaconda.org/bioconda/noarch
                          https://repo.anaconda.com/pkgs/main/linux-64
                          https://repo.anaconda.com/pkgs/main/noarch
                          https://repo.anaconda.com/pkgs/r/linux-64
                          https://repo.anaconda.com/pkgs/r/noarch
          package cache : /opt/conda/pkgs
                          /root/.conda/pkgs
       envs directories : /opt/conda/envs
                          /root/.conda/envs
               platform : linux-64
             user-agent : conda/4.10.3 requests/2.25.1 CPython/3.8.10 Linux/5.10.16.3-microsoft-standard-WSL2 ubuntu/20.04.3 glibc/2.31
                UID:GID : 0:0
             netrc file : None
           offline mode : False

An unexpected error has occurred. Conda has prepared the above report.

```

# Answer
I select the Miniconda installer directly from the official repository:

```
FROM ubuntu:20.04
ENV LANG C.UTF-8
SHELL ["/bin/bash", "-c"]
RUN apt update -qq && apt -qq -y install curl ca-certificates --no-install-recommends
RUN ( cd /tmp/; curl -LO "http://repo.anaconda.com/miniconda/Miniconda3-py38_4.10.3-Linux-x86_64.sh")
RUN bash /tmp/Miniconda3-py38_4.10.3-Linux-x86_64.sh -b -p /opt/conda
ENV PATH=/opt/conda/bin:${PATH}
RUN conda install -y -c conda-forge mamba
RUN mamba create -q -y -c conda-forge -c bioconda -n snakemake snakemake snakemake-minimal --only-deps
RUN echo "source activate snakemake" > ~/.bashrc
ENV PATH /opt/conda/envs/snakemake/bin:${PATH}

```
> 1 votes

# Answer
Try

```
FROM ubuntu:latest
MAINTAINER Conda Development Team <conda@continuum.io>

RUN apt-get -qq update && apt-get -qq -y install curl bzip2 \
    && curl -sSL https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -o /tmp/miniconda.sh \
    && bash /tmp/miniconda.sh -bfp /usr/local \
    && rm -rf /tmp/miniconda.sh \
    && conda install -y python=3 \
    && conda update conda \
    && apt-get -qq -y remove curl bzip2 \
    && apt-get -qq -y autoremove \
    && apt-get autoclean \
    && rm -rf /var/lib/apt/lists/* /var/log/dpkg.log \
    && conda clean --all --yes
RUN conda install -y -c conda-forge mamba
RUN mamba create -q -y -c conda-forge -c bioconda -n snakemake snakemake snakemake-minimal --only-deps
ENV PATH /opt/conda/envs/snakemake/bin:${PATH}

```

The && approach is better as it will minimise the container size (which is going to be big in any case). Then look at an alpine version once its working.

To be honest, I'm not sure about not activating conda base and personally I'd look to correct that.

Alternatively ...

```
docker pull continuumio/miniconda3
docker run -i -t continuumio/miniconda3 /bin/bash

```

Just for a chuckle, a cool approach is

```
docker pull frolvlad/alpine-miniconda3

```

This is super small.

The docker file being

```
FROM frolvlad/alpine-glibc:alpine-3.14

ARG CONDA_VERSION="4.7.12.1"
ARG CONDA_MD5="81c773ff87af5cfac79ab862942ab6b3"
ARG CONDA_DIR="/opt/conda"

ENV PATH="$CONDA_DIR/bin:$PATH"
ENV PYTHONDONTWRITEBYTECODE=1

# Install conda
RUN echo "**** install dev packages ****" && \
    apk add --no-cache --virtual .build-dependencies bash ca-certificates wget && \
    \
    echo "**** get Miniconda ****" && \
    mkdir -p "$CONDA_DIR" && \
    wget "http://repo.continuum.io/miniconda/Miniconda3-${CONDA_VERSION}-Linux-x86_64.sh" -O miniconda.sh && \
    echo "$CONDA_MD5  miniconda.sh" | md5sum -c && \
    \
    echo "**** install Miniconda ****" && \
    bash miniconda.sh -f -b -p "$CONDA_DIR" && \
    echo "export PATH=$CONDA_DIR/bin:\\\$PATH" > /etc/profile.d/conda.sh && \
    \
    echo "**** setup Miniconda ****" && \
    conda update --all --yes && \
    conda config --set auto_update_conda False && \
    \
    echo "**** cleanup ****" && \
    apk del --purge .build-dependencies && \
    rm -f miniconda.sh && \
    conda clean --all --force-pkgs-dirs --yes && \
    find "$CONDA_DIR" -follow -type f \( -iname '*.a' -o -iname '*.pyc' -o -iname '*.js.map' \) -delete && \
    \
    echo "**** finalize ****" && \
    mkdir -p "$CONDA_DIR/locks" && \
    chmod 777 "$CONDA_DIR/locks"

```

Whoever frolvlad is they certainly know their stuff, its a highly optimised image. This will absolutely nail a light-weight conda file. How well conda will work with Snakemake ... it is well worth trying because the alpine image is 5meg (Ubuntu is 100+ Megs) and I forget how big conda is. I need to do this more for sure, because almost certainly what they have done is streamlined conda to the absolute essentials. Get that working your looping the loop.

Why is this potentially super cool? The problem with alpine is its so minimalist getting anything to work is not easy. `frolvlad` in theory solved this bit via conda, because conda is completely self-contained once its running. Thus its doesn't need any system stuff at all, hence alpine and conda in THEORY would be a fantastic solution.

My honest opinion is you are sort of trying to reinvent the wheel (which sometimes that's a good thing). However, I think you'd find the other stuff just works because someone has sat down with yaml and tested everything endlessly.

The problem with your file is you've got to clean it without additional layers being dumped onto it and that's why loads of && because that enables a clean without a layer collapsing onto it blocking the clean. The docker file you create has loads of unnecessary stuff, which you can't delete. I guess your priority is to get it working.
> 0 votes

---
Tags: snakemake, docker
---


title: Trinity de novo assembly not completing

text:
# Question
Title: Trinity de novo assembly not completing
I'm trying to assemble some transcriptomes using Trinity, but am having issues getting Trinity to finish. I've been trying to get this to work for weeks, but am hitting a wall and would very much appreciate any and all help.

I have 2 treatment conditions, with 3 replicates per condition, and they were sequenced on a NovaSeq (using 2x150 base reads). I originally tried assembling all the samples together, but even after running for 8+ days it did not finish.

I am now trying to run just 1 sample from each condition, but even after 4 days Trinity will not finish. It seems like it is stuck at Butterfly, and I'm unsure how to proceed. Some more detail and things I've tried:

These are microbial isolates grown in axenic culture, not metagenomic samples. The microbes are protists though (not bacteria), so their transcriptomes are more complex than prokaryotes, but not as complex as higher eukaryotes like humans/mice.

I do not need to have high sensitivity for lowly expressed transcripts, so I tried setting the "min\_kmer\_cov" higher. First I set it to 2, which did not help. Then I tried 10 (seemed extreme but I was desperate, but even that has not solved it).

-I am running these on a cluster, and admittedly I am not knowledgeable at all about how to optimally request resources, so I am wondering if that is part of the issue. Our cluster uses slurm, and my resource requests right now are

```
#SBATCH --mem 100000
#SBATCH -n 24
#SBATCH -N 6".

```

-I tried running just one replicate from one of my conditions. This **did** finish (after ~42 hours), so it seems like the issue is not an inherent problem with the transcriptomes/the data itself.

Thank you very much!

---

By efficiently constructing and analyzing sets of de Bruijn graphs, Trinity fully reconstructs a large fraction of transcripts, including alternatively spliced isoforms and transcripts from recently duplicated genes. Compared with other de novo transcriptome assemblers, Trinity recovers more full-length transcripts across a broad range of expression levels, with a sensitivity similar to methods that rely on genome alignments. This approach provides a unified solution for transcriptome reconstruction in any sample, especially in the absence of a reference genome.

---

Trinity Trinity is a de novo assembler that fully reconstructs a large fraction of transcripts, including alternatively spliced isoforms and transcripts from recently duplicated genes using de Bruijn graphs. Trinity recover full-length transcripts across a broad range of expression levels, with a sensitivity similar to methods that rely on genome alignments.

Grabherr, MG et al. (2011) Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nature Biotechnology. 2011 May 15;29(7):644-52.

# Answer
> 1 votes
I'm not familiar with the use of Trinity for microbial assembly. There may be some assumptions for animal assembly that are broken in protists, although I can't think of any reason why that would slow the butterfly step down.

It's possible that your dataset is simply too large. One thing you could try is reducing the dataset size (after pooling all samples together). Trinity can do a kmer-based normalization, more details here:

https://github.com/trinityrnaseq/trinityrnaseq/wiki/Trinity-Insilico-Normalization

> Large RNA-Seq data sets, such as those exceeding 300M pairs, are best suited for in silico normalization prior to running Trinity.
>
> Note, all recent versions of Trinity will perform in silico normalization by default. You can turn it off with `Trinity --no_normalize_reads`.
>
> If for any reason, you'd like to normalize data separately, you can do so as follows:
>
> ```
>  %  $TRINITY_HOME/util/insilico_read_normalization.pl
>
> ```

If that normalization is still failing to reduce the dataset down to a reasonable level, it suggests that the dataset is probably more complex than expected. This could be caused by contamination with other microbes, and may be improved by first filtering out contaminant reads.

In any case, if you're having difficulty with getting Trinity to work, the best approach is to contact the developers by logging an issue on the Trinity github:

https://github.com/trinityrnaseq/trinityrnaseq/issues

You could try using an alternative short-read genomic assembler like SPAdes. I see that SPAdes also has an RNA mode and a metaviral mode, so it might already be a bit more closely tuned to smaller transcriptomes:

https://cab.spbu.ru/software/spades/

https://github.com/ablab/spades/blob/spades\_3.15.4/README.md#metapv

---
Tags: rna-seq, trinity
---


title: UPGMA computation

text:
# Question
Title: UPGMA computation
The unweighted pair group method with arithmetic mean (UPGMA) is a hierarchical clustering method, for example used in phylogeny. In a phylogenetic tree it would result in a global molecular clock for taxa under consideration.

Two clusters $A$ and $B$ with the minimum distance are selected and merged to form a new cluster $C$, where the distance between $C$ = $A \cup B$. Assuming another $D$, we can compute: $$ d(C, D) = \frac{1}{|C| \cdot |D|} \sum\_{x \in A \cup B,\\, y \in D} d(x, y). $$

I think $d(x,y)$ is the distance between objects $x$ and $y$ in the input matrix.

Given the information, is there a way to prove/disprove the following statement?

$$ d(C, D) = \frac{1}{|A| + |B|} \Big( |A| d(A, D) + |B| d(B, D) \Big). $$

# Answer
> 0 votes
The key words in UPGMA are 'pairwise' and 'mean'.

In phylogeny the format is a bit unusual because $C$ here is the MRCA (most recent common ancestor) for $A$ and $B$. In other words this represents a three taxa tree. I agree in UPGMA this works, because of its use of the mean (which is also its criticism - below).

Under UPGMA - in your equation $d(C) = \frac{d(A,B)}{N}$ ...where $N=2$, i.e. mean.

Thus under UPGMA

$d(C,D) = d(A,D) - \frac{d(A,B)}{2}$

There is no need to include $d(B,D)$ because $d(B,D) = d(A,D)$, its the mean that was resolved in establishing $C$. If you want to assign $E = D U A$ that easy too.

Your equation doesn't work because |A| and |B| are both 0 under UPGMA, i.e. is pairwise so d(A,A) is always zero. Thus under your equation $d(C,D) = 0$ which simply means $D$, $A$, $B$ are identical and results in a matrix of 0.

---

UPGMA isn't really used in general in biological systems, because the use of the mean rarely makes it a good assumption. I 100% agree has been used in genomics in the modern era, absolutely true, its use is not encouraged because in results in artifactual clusters due its inability to accommodate unequal variation in the underlying genetic process. It could be used in for example in agricultural systems because these conform to uniform distributions, e.g. the distances between apple trees in an orchard. In the natural environment species clump and don't show uniform distributions. UPGMA was used in phylogenetic but the assumption of a global clock very rarely held, so it was discontinued.

---
Tags: phylogenetics, phylogeny, statistics, computation
---


# Question
Title: Workflow for identification of new splicing variants from RNA-seq data?
I'd like to perform a search for possibly unidentified splicing variants of a specific protein in *A. Thaliana*. I do not have my own RNA-seq data. It has been suggested to me to use this workflow:

1. Download RNA-seq data from SRA with `fastq-dump`.
2. Construct an index of reference sequences with `bowtie2-build`.
3. Map the reads from RNA-seq data to the reference sequences with `bowtie2`.
4. Remove non-specific mappings with `samtools view -q 10`. Sort the mappings with `samtools sort`.
5. Browse the aligned reads with `tablet`.

I have been given example data of another organism (not *A. Thaliana*) and have successfully completed these steps. However, I don't know how to apply that in my research. Specifically, I'd like to know:

1. Is there other/better way to do do this?
2. Can I narrow the search/alignment to a specifc gene?
3. Where can I download *A. Thaliana*'s reference genome in FASTA format?
4. Which SRA runs should I choose for my analysis?
5. How to interpret the results to get the information on the different splicing variants?

Linking a similar question I asked before: How to find an unidentified splicing variant of a protein?

# Answer
> 2 votes
> Is there other/better way to do do this?

There is almost always a better way for pretty much anything but in this case I would stick with what your colleagues are using/suggesting. Once you have a worked-out case, you can always improve if needed. Other than that what is suggested to you seems reasonable to me (I have done sth similar for a gene family of some 100 genes using a workflow similar to what was suggested to you)

> Can I narrow the search/alignment to a specifc gene?

Technically yes, instead of the whole genome, you can use a subset of the genome. This will present some problems, for example since the aligner will not be able to find native "targets" for most of the reads, some reads will map to your "new reference" (genome subset) in a "False-positive" manner. You will need to play with the settings. I would advise against this.

> Where can I download A. Thaliana's reference genome in FASTA format?

Genbank, EMBL, ... There might be a specific repo for your organism like there is Flybase for Drosophila species and Xenbase for Xenopus species.

> Which SRA runs should I choose for my analysis?

RNA-seq experiments conducted on your species of interest and ideally on a tissue where your genes of interest are highly expressed. If there is none, you could use a similar organism but this would come with some problems. The most ideal would be RNA-seq with long-read sequencing (like Pacbio).

> How to interpret the results to get the information on the different splicing variants? When you visualize the aligned reads and compare these with existing gene annotation model, it will immediately be obvious if there are novel isoforms, could be due to "new" or not-yet-known exons/splice signals.

---
Tags: rna-seq, rna-splicing
---


title: Mapping to assemblies in subdirectories

text:
# Question
Title: Mapping to assemblies in subdirectories
I am performing a bash loop one liner:.

```
for k in */Assembly/*/*.fastq; do minimap2 -ax map-ont assembly.fasta $k > mapping.sam; echo $k; done

```

The file hierarchy is as follows:

```
\Drain_4
    \Assembly
        \A4_assembly
           assembly.fasta
           A4.fastq
        \D3_assembly
        \I8_assembly
        \U7_assembly
\Drain_5
\Drain_7

```

Where `\` denotes a directory. This is the error i am getting:

```
[ERROR] failed to open file '.fasta': No such file or directory A4.fastq

```

Both input files are in the same directory.

The fastq is my current directory. I'm guessing i will need to provide the complete pathname for both my input files?

I should mention that in all folders the assembly.fasta file is named the same but the fastq file is different (has a different prefix depending on which folder it is in.

# Answer
The most likely issue is that the `assembly.fasta` file is not found in the current directory. Taking @m's comments into consideration, a minimal rewriting of your script to fix this would be as follows:

```
for k in */Assembly/*/*.fastq; do minimap2 -ax map-ont $(dirname ${k})/assembly.fasta $k > $(dirname ${k})/mapping.sam; echo $k; done

```

But I'd recommend at least adding additional quote and brace protection, just in case file names contained spaces:

```
for k in */Assembly/*/*.fastq;
  fDir = $(dirname ${k});
  do minimap2 -ax map-ont "${fDir}/assembly.fasta" "${k}" > "${fDir}/mapping.sam";
  echo ${k};
done

```

edit: fixed up variable capitalisation
> 2 votes

# Answer
The solution is quite simple you issue the one liner either from within `A4_assembly` directory or else `NC_D2` (see below). It will then find `assembly.fasta` and the error will disappear. That's because the current working directory is assigned to the directory in which the one liner is issued. Directory A4\_assembly is where `assembly.fasta` is situated, thus is will work.

If will fail to find the `fastq` file thus you will need to rewrite the loop, simple,

```
for k in *.fastq; do minimap2 -ax map-ont assembly.fasta $k > mapping.sam; echo $k; done

```

... will work, executed from within `Assembly_4`

The one liner suggests you are executing the command from the directory above the `Drain_4` directory. If this is how you intended the one liner to operate I don't think its the solution to the one liner, because there is a risk you will encounter logical errors (not bugs) ..

```
*/Assembly/*/*.fastq

```

If you did execute this commend from above Drain\_4 it will looking all over the place here in every Drain and every assembly directory therein. If you executed this anywhere else it will fail to find \*.fastq because the start of the wildcards is the current working directory in a bash for relative paths.

However, when you dump the file it goes into one one place:

```
 > mapping.sam

```

So there is a risk are performing loads and loads of minimaps but they are simply overwriting each of other, i.e. they are replacing each other.

The other issue is

```
NC_D2/

```

This directory is not represented in your directory file.

Fundamentally the problem is that `assembly.fasta` is always in one location which appears to be in the directory above the directory structure you described, viz. `*/Assembly/*/*.fastq` ... according to the one liner, but in the directory structure its in `A4_assembly` it is also likely to be in `\D3_assembly` amongst loads of others and other Drains (master directories).

Thus even if it 'finds' the reference genome (if that is what it is?), it will see that for everything ... and I mean everything.

There are dynamic solutions to what I think the code using bash is looking to achieve because you can issue regex from within the do ... done, or simply do it as a script.

---

I suggest to simply execute the first script I wrote in every Assembly directory then everything will simply work and there's no risk of logical errors.
> 0 votes

---
Tags: bash, linux, loop
---