# Question
Title: 16S species-level taxonomic assignment--what is the current state of the art?
We currently use DADA2 for picking ASVs and the `assignTaxonomy` funciton for assingment to genera. Google does bring up various recent articles on species-level assignment with 16S.

Does anyone have an opinion on the current best-practice?

More simply, are there any current candidates for reasonably accurate species-level assignment?

# Answer
I have been using dada2 and phyloseq (phangorn for trees) combination to analyze a V4 region 16s rRNA data set. I find their approach and explanations very understandable and usable since I prefer R.

I was under the impression that even full-length 16s rRNA gene sequences were inadequate to classify up to species level since *closely related species within a genus can have identical 16S rRNA gene sequences*.

But, turns out, intra-genomic variants of this gene within a species or even strain may allow their unambiguous identification according to this paper:

Johnson, J.S., Spakowicz, D.J., Hong, B. et al. Evaluation of 16S rRNA gene sequencing for species and strain-level microbiome analysis. Nat Commun 10, 5029 (2019). https://doi.org/10.1038/s41467-019-13036-1) available at: https://rdcu.be/b5cgB

Another relevant article regarding using black box ML algorithms (random forests) versus a straightforward logistic regression is discussed in this one:

Topçuoğlu BD, Lesniak NA, Ruffin MT 4th, Wiens J, Schloss PD. A Framework for Effective Application of Machine Learning to Microbiome-Based Classification Problems. mBio. 2020;11(3):e00434-20. Published 2020 Jun 9. doi:10.1128/mBio.00434-20

Since the OP has probably already moved on long since the original query, I am writing this to document relevant current references for myself and anyone researching this topic on StackExchange Bioinformatics.
> 1 votes

# Answer
My advice is to go for Fiannaca et al. 2018 or similar and steer clear of the DADA2 classification.

My reasoning,

... my criticism of the DADA2 assignments is the reliance on an old machine learning algorithm (naive Bayes). In its favour the taxonomic assignment looks broad. Modern ML will never rely on a single ML method (e.g. lasso regression, SVM, k-means), the only possible exception is "Random Forests"/decision trees, which are popular - but even then. Naïve Bayes (which is not to be confused with 'Bayesian') was popular about a decade ago, which makes sense with the publication cited. The problem is that ML methods give different results.

In its favour it does look well trained, but you would need to understand the distribution used in the training.

---

Again my advice is to go for a deep learning solution, such as Fiannaca et al. 2018 which is trained on bacterial 16s. What looks good here is the CNN architecture. They also implement Deep Belief, which I've never heard of, but suspect it is a variant of RNN. There description suggests they don't much know about it either. Finally they do implement a naïve Bayes solution for comparison.

CNN is the top performing method, which is no surprise. Alot will depend on their filtering and striding, as well as how they trained it. I don't see their training graphs, nor their computer architecture (which are super important). We trust they have done them well.

Ultimately, classification down to the genus level is 85-90%, which for this stuff must be good (for Bioinformatics), but in reality its low.... in fact for DL 85% is not great and needs further attention (when I finally get my CNN working I might live to regret this statement).

Phylogeny (maximum likelihood) is the only error-free way of doing this. There is a bedrock of historic work here. If you want to get down to the species level then trees are your only way to get there. Obviously you'd have a mixed strategy towards that goal.

---

**Final note**, whilst I am critical of DADA2 for classification their variant calling looks very cool, i.e. the robustness of the NGS data. I definitely get the impression the ML is just an add on for them.
> 0 votes

# Answer
DADA2 implements its own method-asignSpecies-for doing exactly this, which uses exact matches to species in a reference database and will be as accurate as any current method - for a single hyper variable region (which means not *hugely* accurate). Basically there's not enough variability in any single hyper variable region to separate species across all genera.

Recent advances on species/strain identification, as far as I'm aware, have been concentrating on using the whole of the 16S region from 3rd generation long read technology.
> 0 votes

---
Tags: fasta, phylogenetics, phylogeny, machine-learning, bacteria
---


# Question
Title: Question: Bismark methylation for non-uniquely mapping BS-seq reads
I am looking for the DNA methylation, specifically inside small RNAs; piRNA. I am mapping BS-seq reads on the genome with bismark. From the user guide of bismark, it maps reads uniquely on the genome, which mean it reports first read with best score and discard rest of the alignment and if a read mapping to multiple location with same score it will discard all such reads.

Typical command for obtaining uniquely mapped reads from bismark

```
~/Software/Bismark-0.22.3/bismark GenomePATH/Original_Genome/  --bowtie2 -p 10 -1 Methylation_reads/Sample1_R1.fq.gz -2 Methylation_reads/Sample1_R2.fq.gz --bam -o Sample1

```

Since, I am interested in piRNA only which may have multiple source in the genome (one piRNA coded from different location of genome), I wanted to take all mapped locations in consideration. Can you suggest what changes I have to make in bismark command. After reading bismark user guide I am confused because I could not able to find relevant parameter.

Also, please let me know if this approach is correct, if not what may be best approach for this.

# Answer
> 1 votes
I'm not certain about your piRNA approach, maybe others will be able to help with that question. As for multi-mapped reads, from Bismark's documentation you can find the `--ambig-bam` flag:

> For reads that have multiple alignments a random alignment is written out to a special file ending in .ambiguous.bam. The alignments are in Bowtie2 format and do not any contain Bismark specific entries such as the methylation call etc. These ambiguous BAM files are intended to be used as coverage estimators for variant callers.

It doesn't look like Bismark will return all of the multiple alignments, just a randomly selected one for a given read or read pair (this probably results from the way Bowtie2 handles multiple alignments within Bismark).

---
Tags: ngs, read-mapping, methylation
---


title: awk field separtor

text:
# Question
Title: awk field separtor
I have a file like this :

and I want the output as:

I tried using awk :

```
awk 'BEGIN { OFS=FS="\t" } { sub("\\..*", "", $1); print }'

```

and it prints the first id ENST. But when I change it to:

```
awk 'BEGIN { OFS=FS="\t" } { sub("\\..*", "", $5); print }'

```

it doesnt work.

What is wrong that I am doing here. Kindly help.

# Answer
> 2 votes
Assuming you want to only show the fifth id from the first field, the following Perl script should work:

```
perl -pe 's/^(.*?)\t//; @F=split("\\|", $1); print(@F[4]."\t")' input.tsv > output.tsv

```

Explanation of the Perl code:

1. Loop over the lines in the input file (this is enabled by the command-line option `-p`)
2. Remove everything up to the first tab character from the input line, and save the removed part (excluding tab character) into the grouping variable `$1` \[identified by brackets\].
3. Split variable `$1` at the pipe character, `|`, saving the split components into the array `@F`
4. Print the fifth field \[0,1,2,3,**4**\] from the array `@F` \[without a line break at the end\]
5. Print the remainder of the input line (this is also enabled by the command-line option `-p`)

# Answer
> 1 votes
You need the first column's 6th component when split by `|` and then the remaining columns of the file where `\t` is the delimiter. Use awk.

```
awk -F "\t" -vOFS="\t" 'NR==1 {print} NR>1{split($1,arr,"|"); print arr[6],$2,$3,$4,$5}' your_file

```

# Answer
> 0 votes
You could just do this in libreoffice calc using text to columns (data menu)

If you want to use awk (or maybe gawk) you could use capture groups; but rather that mess around with a fairly complex regex I'd probably do it in two stages like below

```
awk -F"|" '{print $5,$0}' OFS="\t" your_file|\
awk -F"\t" '{print $1,substr($0,index($0,$3))}'

```

The second awk prints everything except field 2

---
Tags: awk
---


# Question
Title: Determining significance of a variable in a glm model
How would one determine the significance of a variable in a glm model?

If I, for example, have a dataframe like the one seen below, how would I determine if the origin of the sample has a significant effect on the value? (this is the number of enzymes capable of degrading the substrate of that matter)

```
Substrate    variable value origin
cellulose       M09    8    free
mannan          M12    2    free
glycogen        M65    2    free
chitin          M87    4    free
cellulose       M90    2    isolate
manan           M78    1    isolate
glycogen        M21    4    isolate
chitin          M21    1    isolate

```

So far I have tried:

```
mcomp = glm.nb(value ~ origin, data = my_data)

summary(mcomp)
Deviance Residuals:
    Min       1Q   Median       3Q      Max
-0.9625  -0.9047  -0.9047   0.1212   3.5232

Coefficients:
              Estimate Std. Error z value Pr(>|z|)
(Intercept)   -0.01657    0.06571  -0.252  0.80097
originisolate -0.21911    0.08180  -2.679  0.00739 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(0.3418) family taken to be 1)

    Null deviance: 2053.5  on 2679  degrees of freedom
Residual deviance: 2046.3  on 2678  degrees of freedom
AIC: 6517.5

Number of Fisher Scoring iterations: 1

              Theta:  0.3418
          Std. Err.:  0.0186

 2 x log-likelihood:  -6511.4590

```

So free becomes the intercept and then isolate if significantly different from that. Does this mean Origin has a significant effect on the value?

Would the better approach be to do the following?:

```
mcomp = glm.nb(value ~ origin + Substrate, data = comb_data)
summary(aov(mcomp))
              Df Sum Sq Mean Sq F value Pr(>F)
origin         1     23   22.55   6.612 0.0102 *
Substrate     44   1445   32.84   9.631 <2e-16 ***
Residuals   2634   8981    3.41
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’

```

This shows me that origin and substrate have an effect on value if I understand correctly?

# Answer
Second viewing of the question from what I can see -0.22 as a coefficient of origin is a strong negative association, so yeah it has a major impact. Its not how I would have done it, but that looks to be the result.

---

First viewing,

I'm going to throw my hat in here. We don't know what 'origin' is about, anyway just throw everything, i.e. each substrate and the origin into the same regression calculation. Check for a low-residual and preferably do a Q-Q plot, transform your data it this doesn't look good.

The key and the thing you are missing is your regression weights, without that I couldn't say very much. If the regression weight is near zero for 'origin' then it has zero impact. If the regression weight of 'origin' is positively greater than everything else ... I assume there are skewed distributions of 'substrates' between the 'origins'. If the regression weight of 'origin' is negative but still greater than all other regression weights then it is adversely affecting the 'value' you are seeking.

I don't know the experiment, the biological system or really the 'substrate' assays, so I can't comment any further.

The two issues I have are:

1. Doing an ANOVA on the output of a regression analysis doesn't make much sense to me. It is not something I would do, nor something in ML or GLM I've encountered.
2. Are you doing pairwise substrate/origin calculations? I presue not, but just in case this not how GLM works.
> 1 votes

# Answer
There is no better method, it's a matter of what you want to test or what is your question.

Using the `anova()` or `aov()`, test the terms collectively. For example, in your example with Substrate, the null hypothesis is that the coefficients are all zero, meaning `cellulose =0, mannan =0 , ....`

If the question is, "do the isolate samples have a higher value than origin samples?", then you can use your first model, where `free` is set as the reference and you test whether the effect of `isolate` is non-zero. Likewise you can do this for substrate and set of them as your reference. You can also do other pairwise comparisons using this model.

If the question is, "does origin have a significant effect on value, after controlling for substrate?", then you can use your second model.
> 0 votes

---
Tags: modelling, linear-regression, statsmodels, log-likelihood
---


title: which is the right prosite pattern?

text:
# Question
Title: which is the right prosite pattern?
I would like to ask if the right Prosite pattern for this multiple alignment:

Is it this:

```
A-T-[AT]-G-x-C-[AGC]-C-x(1,4)-A

```

or this:

```
A-T-[AT]-G-x-C-[AGC]-C-x(1,3)-[ATGC]-A

```

Actually, I would like to ask if `x` always implies a gap, and if it is proper to use it when any symbol but not gap can be in the sequence. Thank you!

# Answer
> 1 votes
If nucleotides are acceptable in prosite format because I thought it was protein motifs, the most accurate would be,

```
A-T-[AT]-G-x-C-[AGC]-C-[GCA](0,1,3)-x-A

```

The key point is that,

1. `[GCA](0,1,3)` of `x(0, 1, 3)` means zero/ nucleotides (3 indels) or one/ G or C or A nucleotide, and 2 indels or three/ G or C or A nucleotides.
2. `x` is any nucleotide i.e. A or G or T or C - it doesn't mean an indel. `x(0, 1)` is the best representation of a single indel (you call them gaps), here meaning any nucleotide may or may not be present at a given site.

Example from Prosite documentation:

> `<A-x-[ST](2)-x(0,1)-V`
>
> This pattern, which must be in the N-terminal of the sequence ('\<'), is translated as:
>
> Ala-any-\[Ser or Thr\]-\[Ser or Thr\]-(any or none)-Val

---
Tags: sequence-alignment, prosite
---


title: Modify how a script prints the output?

text:
# Question
Title: Modify how a script prints the output?
I'm trying to print out all measured distances between a ligand and 2 residues (using PyMol), below is the script that I'm using. I only want to measure the distances between all the ligand's O & ND2 of N312, as well as all the ligand's N & OD2 of D113.

```
from pymol import stored
stored.a, stored.b, stored.c, stored.d = [],[],[],[]
cmd.iterate_state(1, "r55_01///UNL`1/N*", "stored.a.append(ID)")
cmd.iterate_state(1, "r55_01///UNL`1/O*", "stored.b.append(ID)")
cmd.iterate_state(1, "r55_01//R/ASP`113/OD2", "stored.c.append(ID)")
cmd.iterate_state(1, "/r55_01//R/ASN`312/ND2", "stored.d.append(ID)")

f=open("dst4.txt", 'w')

for a in stored.a:
  for b in stored.b:
    for c in stored.c:
      for d in stored.d:
        f.write( "distance from %s to %s is %s\n" % (a, c,
cmd.distance( "id %s" % a, "id %s" % c,cutoff=5,mode=2)))
        f.write( "distance from %s to %s is %s\n" % (b, d,
cmd.distance( "id %s" % b, "id %s" % d,cutoff=5,mode=2)))

f.close()

```

The output file currently looks like this

```
distance from 9 to 1334 is 4.589359760284424
distance from 8 to 4204 is 3.1509475708007812
distance from 9 to 1334 is 4.589359760284424
distance from 15 to 4204 is 0.0
distance from 9 to 1334 is 4.589359760284424
distance from 17 to 4204 is 0.0
distance from 9 to 1334 is 4.589359760284424
distance from 20 to 4204 is 0.0
distance from 23 to 1334 is 0.0
distance from 8 to 4204 is 3.1509475708007812
distance from 23 to 1334 is 0.0
distance from 15 to 4204 is 0.0
distance from 23 to 1334 is 0.0
distance from 17 to 4204 is 0.0
distance from 23 to 1334 is 0.0
distance from 20 to 4204 is 0.0

```

I want my output to look more like this. Thanks in advance!

```
Lig_N D113 4.59
Lig_O N312 3.15

```

# Answer
> 1 votes
Let's simplify...

`cmd.distance(xx,xxx)` accepts two selector strings and is fine with `and state 1`. And best of all if you select multiple atoms it gives you the smallest distance!

So all you need to do is:

```
# measure
d_O = cmd.distance('UNL and element O', 'resi 312 and name ND2')
d_N = cmd.distance('UNL and element N', 'resi 113 and name OD2')
# write
f=open("dst4.txt", 'w')
f.write(f'Lig_N D113 {d_N:.2f}\n')
f.write(f'Lig_O N312 {d_O:.2f}\n')
f.close()

```

If you are using python 2 (if so do note pymol scripting is better in 3), then the second bit is

```
f=open("dst4.txt", 'w')
f.write('Lig_N D113 %.2f\n' % (d_N,))
f.write(f'Lig_O N312 %.2f\n' % (d_O,))
f.close()

```

## Unsolicited coding tips

A few unsolicited tips... Your original script is rather over-complicated and obfuscated and uncommented, so for better legibility why not do:

* `cmd.iterate_state` is always harder to follow than the `for atom in get_model('xxxx').atoms` synthax
* unless you are stuck with python 2, use f-strings not the `%s` sprintf style formatting
* the `/x/x/x/x` selector synthax is never as readable as the longer `resi`/`resn`/`name` selection algebra.

Code golf is always tempting, but it can result in future-you getting stuck as happened here...

---
Tags: python, pymol
---


# Question
Title: Is it possible to predict protein-ligand binding kinetics using machine learning?
I would like to work on a project that involves the prediction of protein-ligand binding kinetics. What might be the feature that is relevant for the prediction?

# Answer
> 1 votes
> This is a comment to redirect to other answers, but go too long.

Docking is the *in silico* prediction of where a ligand binds in a protein.

There are several Q&As about docking here that are relevant, so please check out:

---
Tags: proteins, machine-learning, interactions
---


# Question
Title: How can I obtain a list of all NCBI gene ID's along with their full name, symbol, and also known as symbols?
I would like to do two things:

(1) Obtain a list of all NCBI (Entrez) gene ID's.

(2) Obtain the "Official Symbol", "Official Full Name", and "Also known as" fields from the NCBI gene pages for each of those genes.

This information, since it is displayed on the NCBI website, must be located in some databases, but I have not had luck finding a tool that will allow me to extract this information.

# Answer
> 3 votes
All of these data are in the `gene_info.gz` file on the NCBI FTP site: https://ftp.ncbi.nlm.nih.gov/gene/DATA/ Specifically, you may want to look at the `gene_info.gz` file. It has the following fields (shown for an example gene below):

```
                              #tax_id [  1]: 9606
                               GeneID [  2]: 5768
                               Symbol [  3]: QSOX1
                             LocusTag [  4]: -
                             Synonyms [  5]: Q6|QSCN6
                              dbXrefs [  6]: MIM:603120|HGNC:HGNC:9756|Ensembl:ENSG00000116260
                           chromosome [  7]: 1
                         map_location [  8]: 1q25.2
                          description [  9]: quiescin sulfhydryl oxidase 1
                         type_of_gene [ 10]: protein-coding
   Symbol_from_nomenclature_authority [ 11]: QSOX1
Full_name_from_nomenclature_authority [ 12]: quiescin sulfhydryl oxidase 1
                  Nomenclature_status [ 13]: O
                   Other_designations [ 14]: sulfhydryl oxidase 1|quiescin Q6 sulfhydryl oxidase 1|testis tissue sperm-binding protein Li 62n|thiol oxidase 1
                    Modification_date [ 15]: 20200604
                         Feature_type [ 16]: -

```

If you are starting with a set of NCBI GeneIDs and want this sort of information for just those, you can use Entrez Direct for that. If you are starting with RefSeq transcript or protein accessions, be sure to check out the other files in the FTP path mentioned above.

# Answer
> 1 votes
Go to Ensembl (https://www.ensembl.org/index.html), click `BioMart` on the menu at the top of the page. Select the dataset you want to use, then go to Attributes, select the feature of each gene you want to get. Then click "Results" button to make sure that the table looks like what you want. Then output the results to a file.

---
Tags: gene, phylogenetics, text-processing
---


# Question
Title: R not being able to find the file despite being in correct directory
I have question regarding a project I have to do about a health student study. I cannot import the csv file via read.csv or read.csv2 despite having transformed all files to csv and being in the correct folder. Here are the command lines:

```
setwd ("C:/Users/fek/Desktop/Analyse i-Share")
getwd()
TA<-read.csv2(file=Biobanque.csv, header=TRUE, sep=",", dec=",", fill=TRUE)

```

And here are the results in the console window:

```
setwd ("C:/Users/fek/Desktop/Analyse i-Share")

getwd()

```

\[1\] "C:/Users/fek/Desktop/Analyse i-Share"

```
TA<-read.csv2(file=Biobanque.csv, header=TRUE, sep=",", dec=",", fill=TRUE)

```

\[1\] Error in read.table(file = file, header = header, sep = sep, quote = quote, : objet 'Biobanque.csv' introuvable

To translate from french, it says it can't find the file!

# Answer
> 4 votes
The problem is just what it says. By omitting the quotes around "Biobanque.csv", R thinks you are talking about an **object whose name is Biobanque.csv**. And R is confused because of course you have not created an object by that name.

---
Tags: r
---


# Question
Title: Do I need to fill in missing side chains and loops before protein-ligand docking for VLS?
I am going to screen a set of small molecules against a **CYP450 2C9** protein target. The structure I am using exactly is 5A5I. Do I need to fill in missing side chains and loops before performing VLS? I assume while preparing the protein I should leave water molecules in the structure as there is a bond between $Fe^{2+}$ and water that is involved in the binding of a ligand.

# Answer
> 2 votes
# In general

Most if not all docking applications will tolerate a jump. Furthermore, you really don't want to dock anything to the flexible regions even if you are using a program that does properly protein flexibility. Furthermore, if your protein adopts a closed conformation upon native substrate binding (not really the case here), you should start with that for better docking results.

# Cytochrome P450

So most Cytochrome P450s are strange though and in many those loops play a large part in their catalysis. I am on a paper under review (bioarxiv) that deals exactly with this in a different cytochrome (P450-BM3, the biotech one). Namely even if the active site is fairly large for an active site in the core of an enzyme, the ligand (testosterone) did not really rotate around and instead the specificity of oxidation was controlled by how the ligand got in. So those crystallographically unresolved loops act as gatekeepers for what can bind. In terms of VLS, one aims for blocking off the cofactor, i.e. the heme. So you can do docking without worrying about them in a first instance, but there is a possibility that a lot of active drugs instead interact with the loops. Certainly, dynamic undocking and other MD methods may be required a few iterations in, to see why some candidates failed and why others worked.

An interesting unofficial further possible spin is the fact that solvent-exposed sequential residues Asp-Pro can hydrolyse under acidic conditions, in a mechanism similar to isoaspartate formation, and several cytochromes have these on the entrance loops and have weird kinetics and smeary SDS-PAGE gels.

Lastly, you are totally correct in your observation. The water binding the heme in the crystal is very special. It is taking the place of oxygen. That could be modelled as is, changed to a pre-encounter O2, or event changing the heme to a heme bound to a radical oxygen —note that some forcefields (CHARMM, Rosetta Ref2015) do not like this weird valency iron.

---
Tags: protein-structure, docking
---


# Question
Title: Difference between isoforms and paralogs in transcriptoms?
I've assembled RNAseq data into a transcriptoms using Trinity. There's a option to keep only the longest isoforms for each transcripts and it lead me to wonder how it deals with duplicated genes (--\>transcripts). I haven't found an anwser in the original paper. Althought I reckon paralogs and isoforms are two distinct features I don't see a way to make a distinction when you only have trasncripts.

# Answer
> 1 votes
It depends. If you are performing Illumina based RNAseq, which is likely, it is definitely tricky, You should be able to do this in an ideal world because,

1. a paralogue will contain SNPs, i.e. there will be nucleotide variation between paralogues, but the mRNA can of an identical length
2. isoform will also have identical sequence but the length of the mRNA will vary depending on what initiation codon it used

I sort of agree however because,

1. there will be 'grey-zones' where paralogues are different sizes and VCF/ SNP calling is not perfect and ...
2. assembly of paralogues can be difficult, risking collapsing a paralogue into a single mRNA, but it is not impossible depending on the level of SNP variation
3. assembly of isoforms would be really hard

I personally wouldn't know how point 3 is solvable, but I guess it is if your isoform is relatively small, i.e. well within the read-length of Illumina output.

# Answer
> 3 votes
If a reference genome assembly exists, then reference-guided assembly with Trinity will keep paralogs, because it treats each genomic region as a different assembly set.

If a reference genome doesn't exist, paralogs will probably result in the construction of non-representative "isoforms" containing variants present in multiple genes. Inchworm does a greedy assembly of reads, so if a chimeric transcript *can* be constructed from the reads (which is likely with paralogous genes), then it might end up being fed to the subsequent steps. Whether that is actually emitted at the end by Butterfly will mostly depend on how the paths through the putative assemblies are traversed, which probably depends on the tolerance of Bowtie to errors.

If paralogous transcripts are sufficiently similar that the differences are within the error tolerance for assembly, then I would expect that they'd be collapsed into a single transcript.

---
Tags: rna-seq, transcriptome, isoform, paralogs
---


# Question
Title: how to extract a different substring for each row of a dataframe in R?
I have a dataframe with putative alignments of thousands of probe sequences. I would like to pull out the adjacent nucleotide from the genome for each alignment. My dataframe includes `seq.name` for the target chromosome and I've added an `adjacent_pos` variable with the location on the target chromosome that I want. Using `read.fasta` from the phylotools package gives me a genome in the format of a dataframe with `seq.name` and `seq.text` for each chromosome and unplaced scaffold.

Working on a test data set with a single chromosome, I successfully used something like this:

```
adj <- cbind(alignments$adjacent_pos, alignments$adjacent_pos)
alignments$adjacent_base <- str_sub(genome$seq.text, adj)

```

Now back to the whole genome, I can't wrap my head around how to direct `str_sub` to the correct `seq.name` for each `seq.position`. I'm sure there's an elegant solution, but I haven't been able to find it!

# Answer
Say your file looks like this:

```
cat test.fa
>chr01
ACCGGTT
>chr02
ATGCATGC

```

Something like what you have done, to read in the fasta and extract the sequence. Key part is creating a GenomicRanges object:

```
library(GenomicRanges)
library(BSgenome)
library(Biostrings)

adj <- data.frame(chr=c("chr01","chr02"),start=c(2,3),end=c(2,3))

    chr start end
1 chr01     2   2
2 chr02     3   3

adj = makeGRangesFromDataFrame(adj)

fasta <- readDNAStringSet("test.fa")
SEQ = getSeq(fasta,makeGRangesFromDataFrame(adj))

as.character(SEQ)
[1] "C" "G"

```

or if the fasta file is huge, you can use Rsamtools, indexing the fasta file first:

```
library(Rsamtools)
fasta = FaFile("test.fa")
indexFa('test.fa')

as.character(getSeq(FaFile("test.fa"),adj))
chr01 chr02
  "C"   "G"

```
> 0 votes

---
Tags: r
---


title: Error after trimming illumina adapters

text:
# Question
Title: Error after trimming illumina adapters
I am removing illumina adapters of the NGS data with a loop. My NGS data is storage in /data/HTS\_seq/. I used this function:

```
 for infile in /data/HTS_seq/*_read1.fastq; do outfile=$infile\_trim.fastq cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -e 0.1 -O 3 -m 30 --max-n 2 -o $outfile $infile; done

 for infile in /data/HTS_seq/*_read2.fastq; do outfile=$infile\_trim.fastq cutadapt -a AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA -e 0.1 -O 3 -m 30 --max-n 2 -o $outfile $infile; done

```

and It gives to me this error:

```
  IOError: [Errno 20] Not a directory: '/data/HTS_seq/PV589_03664AAB_AGGACT_read1.fastq/_trim.fastq'
  Traceback (most recent call last):
  File "/usr/local/bin/cutadapt", line 10, in <module>
cutadapt.main()
  File "/usr/local/lib/python2.7/dist-packages/cutadapt/scripts/cutadapt.py", line 652, in main
options.discard_untrimmed)
  File "/usr/local/lib/python2.7/dist-packages/cutadapt/scripts/cutadapt.py", line 327, in trimmed_and_untrimmed_files
trimmed = untrimmed = xopen(output_path, 'w')
  File "/usr/local/lib/python2.7/dist-packages/cutadapt/xopen.py", line 178, in xopen
return open(filename, mode)
      IOError: [Errno 20] Not a directory: '/data/HTS_seq/PV589_03664AAB_AGGACT_read1.fastq/_trim.fastq'

```

Which is my error? Thank you in advance.

# Answer
> 1 votes
First of all, you aren't actually setting the `$outfile` variable. Your loop has this basic format:

```
for var in *; do newVar="something" command $newVar; done

```

For example:

```
$ for var in 1 2 3 4; do newVar="new:$var" echo "$var: $newVar"; done
1:
2:
3:
4:

```

As you can see, `newVar` is never actually being set. That's because setting a variable is its own operation. So you need to do that, end it with a control character like `;`, or `&`, or a newline, and *then* run your `cutadapt`:

```
$ for var in 1 2 3 4; do newVar="new:$var"; echo "$var: $newVar"; done
1: new:1
2: new:2
3: new:3
4: new:4

```

Or, with your data:

```
for infile in /data/HTS_seq/*_read1.fastq; do
    outfile="$infile"_trim.fastq
    cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \
      -e 0.1 -O 3 -m 30 --max-n 2 -o "$outfile" "$infile"
done

for infile in /data/HTS_seq/*_read2.fastq; do
    outfile="$infile"_trim.fastq
    cutadapt -a AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA \
      -e 0.1 -O 3 -m 30 --max-n 2 -o $outfile $infile
done

```

But that would give a weird file name, like `foo_S1_L001_read1.fastq_trim.fastq`. You probably want names like `foo_S1_L001_read1_trim.fastq` instead, so you would want:

```
for infile in /data/HTS_seq/*_read1.fastq; do
    outfile="${infile%.fastq}"_trim.fastq
    cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \
      -e 0.1 -O 3 -m 30 --max-n 2 -o "$outfile" "$infile"
done

for infile in /data/HTS_seq/*_read2.fastq; do
    outfile="${infile%.fastq}"_trim.fastq
    cutadapt -a AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA \
      -e 0.1 -O 3 -m 30 --max-n 2 -o $outfile $infile
done

```

---

Now, the specific error you show makes me think that isn't the command you ran. Could you perhaps have run this one instead by mistake (note the `/` instead of `\` before the `_trim.fastq`)?

```
for infile in /data/HTS_seq/*_read1.fastq; do
 outfile=$infile/_trim.fastq cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \
-e 0.1 -O 3 -m 30 --max-n 2 -o $outfile $infile
done

```

Your error message is:

```
  IOError: [Errno 20] Not a directory: '/data/HTS_seq/PV589_03664AAB_AGGACT_read1.fastq/_trim.fastq'

```

Which means that `cutadapt` tried to open `/data/HTS_seq/PV589_03664AAB_AGGACT_read1.fastq/` as a directory and then look for a file named `_trim.fastq` inside it. The only way I can think that this would happen is if you ran the command I gave above and not the one you show in your question.

# Answer
> -1 votes
There is a bug here

```
/data/HTS_seq/*_read1.fastq

```

This will go to your root directory `/` and there is definiately no folder called `data` in the root. To prove it type `ls /` and that is where the command is pointing, then compare this to `ls ./`.

To solve it go to the directory `HTS_seq` and type `pwd` then copy and paste the directory structure. I suspect `./data/HTS_seq` will solve it

```
 for infile in ./data/HTS_seq; ...

```

This will look for the folder `data` in your current working directory rather than the root directory.

---

Okay you've been assigned /data as root ... no biggie.

```
PV589_03664AAB_AGGACT_read1.fastq/_trim.fastq

```

Why is `/` in the file name? What happens if you `ls /data/HTS_seq/PV589_03664AAB_AGGACT_read1.fastq/_trim.fastq`? and `ls /data/HTS_seq/` Finally could you please `ls /data/HTS_seq/PV589_03664AAB_AGGACT_read1.fastq\/_trim.fastq`?

The `\/` is important. What I think has happened is a `/` is in the file name and that will cause problems because it will think that `_trim.fastq` is a directory when its part of a file name, and would explain the error.

Then do the same thing using `head` instead of `ls` see what happens.

If the last `ls` and `head` works then rename the files by typing `rename` in the command line. If `rename` works it will be easy to remove `/` from all your files using the power of the `rename` function (its a Perl function BTW).

---
Tags: ngs, cutadapt, loop
---


title: Multiple input & output files?

text:
# Question
Title: Multiple input & output files?
I'm trying to print out the measured distances for each .pdbqt file, and then write a new output file for each .pdbqt, using the original name of the .pdbqt file as the new .txt file name. Below is the script.

I'm using python3 & I'm trying to do this via Linux CLI.

```
from pymol import stored
stored.a, stored.b, stored.c, stored.d = [],[],[],[]
cmd.iterate_state(1, "r55_01///UNL`1/N*", "stored.a.append(ID)")
cmd.iterate_state(1, "r55_01///UNL`1/O*", "stored.b.append(ID)")
cmd.iterate_state(1, "r55_01//R/ASP`113/OD2", "stored.c.append(ID)")
cmd.iterate_state(1, "r55_01//R/ASN`312/ND2", "stored.d.append(ID)")

f=open("tt.txt", 'w')

for a in stored.a:
  for c in stored.c:
      d1=cmd.distance("id %s" % a, "id %s" % c,cutoff=5)
      f.write("%s D113 %s\n" % (a,d1))
      print(d1)

for b in stored.b:
  for d in stored.d:
      d2 = cmd.distance("id %s" % b, "id %s" % d, cutoff=5)
      f.write("%s N312 %s\n" % (b,d2))
      print(d2)

f.close()

```

My folder looks something like this

```
.
├── dst6.txt
├── sp
│   ├── r55_01.pdbqt
│   ├── r55_02.pdbqt
│   ├── r55_03.pdbqt
│   ├── r55_04.pdbqt
│   ├── r55_05.pdbqt
│   ├── r55_06.pdbqt
│   ├── r55_07.pdbqt
│   ├── r55_08.pdbqt
│   ├── r55_09.pdbqt
│   ├── r55_10.pdbqt
│   ├── r56_01.pdbqt
│   ├── r56_02.pdbqt
│   ├── r56_03.pdbqt
│   ├── r56_04.pdbqt
│   ├── r56_05.pdbqt
│   ├── r56_06.pdbqt
│   ├── r56_07.pdbqt
│   ├── r56_08.pdbqt
│   ├── r56_09.pdbqt
│   └── r56_10.pdbqt
└── t6.py

```

THANKS IN ADVANCE!!

# Answer
You made the hardest work, you were not far to achieve what you wanted ;)

This script.py should work after you made it executable (in a terminal type: `chmod u+x script.py`), just type ./script.py in your terminal in your working directory (so here you should be in your "sp" directory). Note that I got inspired by: https://www.biostars.org/p/42474/ to know how to launch pymol within the script and also that I would personally not use pymol for doing this kind of analysis/data extraction from a pdb/pdbqt file but rather use a dedicated analysis package (there are plenty like MDanalysis). I personally use pymol (pymol script) for visualization purpose.

```
#!/usr/bin/env python

import __main__
__main__.pymol_argv = ['pymol', '-qc']  # Quiet and no GUI

from pymol import stored
import pymol, os

pymol.finish_launching()

#put what you've done into a function
def myfunc(file):
    spath = os.path.abspath(file)
    sname = spath.split('/')[-1].split('.')[0]

    pymol.cmd.load(spath, sname)

    stored.a, stored.b, stored.c, stored.d = [], [], [], []
    pymol.cmd.iterate_state(1, "%s///UNL`1/N*" % sname, "stored.a.append(ID)")
    pymol.cmd.iterate_state(1, "%s///UNL`1/O*" % sname, "stored.b.append(ID)")
    pymol.cmd.iterate_state(1, "%s//R/ASP`113/OD2" % sname, "stored.c.append(ID)")
    pymol.cmd.iterate_state(1, "%s//R/ASN`312/ND2" % sname, "stored.d.append(ID)")

    f = open('%s.txt' % sname, 'w')

    for a in stored.a:
        for c in stored.c:
            d1 = pymol.cmd.distance("id %s" % a, "id %s" % c, cutoff=5)
            f.write("%s D113 %s\n" % (a, d1))
            print(d1)

    for b in stored.b:
        for d in stored.d:
            d2 = pymol.cmd.distance("id %s" % b, "id %s" % d, cutoff=5)
            f.write("%s N312 %s\n" % (b, d2))
            print(d2)

    f.close()

mypdbqt = os.popen('ls *.pdbqt').read().split()

for i in mypdbqt:
    myfunc(i)

pymol.cmd.quit()

```
> 1 votes

---
Tags: python, pymol
---


# Question
Title: Most scalable pseudotime ordering algorithm
What algorithms for linear Pseudotime trajectory construction (diffusion-based) are the most scalable to large datasets?

I'm currently using Slingshot based on the recommendation in this manuscript: https://www.biorxiv.org/content/biorxiv/early/2018/03/05/276907.full.pdf

Slingshot is nice and easy to use, but it turns out it is terrible at scaling to big datasets.

For 2000 highest variable genes and a random sampling of a variable number of cells in my ~80% sparse expression matrix I get the following computation times in R (96 GB RAM):

```
#cells / Time (min)
100 / 0.7
1000 / 4.3
2000 / 12.3
5000 / 71.1

```

At this rate it would take nearly 1.5 days to process 200,000 cells and 4.5 months to process 2 million cells based on 2,000 HVG.

I feel like Slingshot is limited by a fundamentally unnecessary comparison of diffusion for all cells against all cells, rather than preclustering and then finer resolution of pseudotime across cluster edges.

What algorithms are better, but still implement a similar realization of pseudotime ordering? I just need a single linear trajectory.

# Answer
> 1 votes
As it turns out, pseudotime algorithms are exceedingly difficult to parallelize. Pseudotime does not scale linearly to a subset of data compared to the full dataset, and as such a course-grained approach to parallelizing pseudotime calculations on random bins of the data, followed by bin merging, does not work. These values will never converge on the actual full-dataset solution.

I put in a pretty good effort trying to parallelize slingshot, but with no success. I realize my code may be of limited utility without more explanation, but I'm not putting more time into this. FWIW:

This code:

1. Finds computationally optimal bin size for slingshot based on a number of genes and features
2. Randomly splits the SingleCellExperiment object into that number of bins and runs slingshot on each bin
3. Orients each bin so that a cell with a pseudotime of 1 in the first bin is similar to a pseudotime of 1 in every bin. Pseudotime can get flipped depending on the sampled subset.
4. Merges all bins into a single matrix so each cell has a assigned pseudotime based on randomly sampled bins.
5. Repeats steps #2 to #5 with a different random seed for binning every time.
6. Measures convergence towards actual solution at three pseudotime values, where convergence is defined as correlation of cell indices along the trajectory.

I abandoned the project after observing that more iterations to not change the convergence towards the actual solution. While pseudotime values close to 0 recapitulated the actual solution quite well, pseudotime values towards 100 failed to correlate at all with the actual solution. There was no convergence with more iterations.

Maybe this will be helpful to somebody, which is why I'm posting.

```
# Parallel Iterative Sampling with Slingshot
# pissshot is a course-grained estimator of the actual slingshot solution, and iteratively converges towards the actual solution
# pissshot randomly partitions cells in a large database into bins, bin size is determined by finding the computationally optimal number of cells
# pissshot is run independently on each bin, pseudotime vectors are aligned across all bins, and then all bins are merged
# pissshot repeats this process and the average of pseudotime assignments for each cell is taken at the end of each iteration.
# The convergence of the blingshot model at the end of each iteration is calculated by running slingshot on a small sample of adjacent cells at a random pseudotime interval and measuring correlation with the predicted model
# When satisfactory convergence is reached, blingshot returns a singlecellobject, just as slingshot would

pissshot <- function(sce, num_iterations = 20){
    cat("Running blingshot on", dim(sce)[2],"cells and",dim(sce)[1],"features\n")
    cat("  Step 1/5: Finding computationally optimal increment sizes\n   ")
#    increment_size <- FindOptimalIncrementSize(sce)
    increment_size <- 688
    num_increments <- ceiling(dim(sce)[2]/increment_size)-1
    cat("\n  ...Increment size of",increment_size,"cells is computationally optimal\n")
    cat("  Step 2/5: Course-grained pseudotime assignment across",num_increments,"increments")
    ps <- list()
    convergenceArr <- c()
    for(iteration in seq(from=1,to=num_iterations,by=1)){
        cat("\n   ... iteration",iteration,"...\n")
        cat("         randomly assigning cells to",num_increments,"bins\n         0%.")
        bins <- PartitionBins(sce, increment_size = increment_size, seed = 7*iteration)
        cat("100%\n         running slingshot on each bin\n         0%")
        cell_pseudotimes <- c()
        for(b in seq(from=1,to=length(bins),by=1)){
            bin_with_pseudotime <- suppressMessages(slingshot(bins[[b]]))
            # Find right cell and left cell in pseudotime array (cell with highest/lowest pseudotime values)
            pcells <- as.data.frame(colData(bin_with_pseudotime)$slingPseudotime_1)
            rownames(pcells) <- colnames(bin_with_pseudotime)
            pcells <- data.frame(lapply(pcells, function(x) as.numeric(as.character(x))), check.names = F, row.names = rownames(pcells))
            pcells.sorted <- pcells[order(-pcells[,1]), , drop = FALSE]
            right_cell_ID <- rownames(pcells.sorted)[dim(pcells.sorted)[1]]
            left_cell_ID <- rownames(pcells.sorted)[1]
            right_cell_pos <- match(right_cell_ID,colnames(bins[[b]]))
            left_cell_pos <- match(left_cell_ID,colnames(bins[[b]]))
            right_cell_logcounts <- assays(bins[[b]])$logcounts[,right_cell_pos]
            left_cell_logcounts <- assays(bins[[b]])$logcounts[,left_cell_pos]
            if(b==1&&iteration==1){
                ref_right_cell_logcounts <- right_cell_logcounts
                ref_left_cell_logcounts <- left_cell_logcounts
            }
            # Figure out whether this pseudotime in this bin has flipped relative to reference (first bin), and if so, flip pseudotime values
            cor_opposite_ends <- mean(cor(right_cell_logcounts,ref_left_cell_logcounts),cor(left_cell_logcounts,ref_right_cell_logcounts))
            cor_same_ends <- mean(cor(right_cell_logcounts,ref_right_cell_logcounts),cor(left_cell_logcounts,ref_left_cell_logcounts))
            flipped <- FALSE
            if(cor_opposite_ends > cor_same_ends){
                flipped <- TRUE
            }
            if(flipped==TRUE){
                # find max pseudotime in colData(bin_with_pseudotime)$slingPseudotime
                # Recalculate pseudotime in bin_with_pseudotime as max-value
                pseudotimes <- colData(bin_with_pseudotime)$slingPseudotime_1
                colData(bin_with_pseudotime)$slingPseudotime_1 <- max(pseudotimes)-pseudotimes
            }
            cells_with_pseudotime <- as.data.frame(colData(bin_with_pseudotime)$slingPseudotime_1,colnames(bin_with_pseudotime))
            cell_pseudotimes <- rbind(cell_pseudotimes, cells_with_pseudotime)
            cat(".")
        }
        cat("100%\n")

        ps[[iteration]] <- cell_pseudotimes[order(row.names(cell_pseudotimes)), , drop = FALSE]

        pss <- data.frame(ps)
        if(iteration>1){
            convergence10 <- measureConvergence(pss, sce = sce, ptime = 10, increment_size = increment_size)
            convergence50 <- measureConvergence(pss, sce = sce, ptime = 50, increment_size = increment_size)
            convergence99 <- measureConvergence(pss, sce = sce, ptime = 99, increment_size = increment_size)
            cat("   convergence at 10:",convergence10,"... at 50:",convergence50,"...at 100:",convergence99,"\n")
            convergenceArr <- rbind(convergenceArr,c(iteration,convergence10,convergence50,convergence99))
        }
        # Iterate until the model converges towards the actual slingshot solution, within an indicated fraction (i.e. convergence = 0.02)
        # The convergence measure compares how similar the model is to the actual slingshot solution, i.e. a convergence of 0.02 means that cell pseudotime values are within 2% of the actual solution
        # Assess this by selecting a pseudotime at random, analyzing the following increment_size cells, and comparing accuracy of pseudotime measurement across that trajectory with the composite of the course-grained analysis

        # need to study how convergence changes over pseudotime
    }

    # remove cells with really high standard deviations

    return(pss)
#    return(convergenceArr)
}

measureConvergence <- function(mat, sce = sce, ptime = 45, increment_size = 500){
    # this line is throwing an error
    mat[,"avg"] <- apply(mat[,1:dim(mat)[2]],1,mean)
    mat <- mat[order(mat[,"avg"]), , drop = FALSE]
    # select a random pseudotime increment
    mat <- mat[mat[, "avg"] > ptime,]
    mat <- mat[1:increment_size,]
    cell_IDs <- rownames(mat)
    # pull out this list of cell_IDs from the sce object and run slingshot
    sce2 <- suppressMessages(slingshot(sce[,cell_IDs]))
    cpcells <- as.data.frame(colData(sce2)$slingPseudotime_1)
    rownames(cpcells) <- colnames(sce2)
    cpcells <- data.frame(lapply(cpcells, function(x) as.numeric(as.character(x))), check.names = F, row.names = rownames(cpcells))
    # sort cpcells by increasing pseudotime
    cpcells <- cpcells[order(cpcells[,1]), , drop = FALSE]
    ordered_cells <- cbind(rownames(mat),rownames(cpcells))
    # get the index of each cell in the dataframe, run a correlation on how well the indices line up
    cell_indices <- c()
    for(cell in cell_IDs){
        new_indices <- c(which(ordered_cells[,1] == cell), which(ordered_cells[,2] == cell))
        cell_indices <- rbind(cell_indices,new_indices)
    }
    convergence_measure <- cor(cell_indices[,1],cell_indices[,2])
    plot(cell_indices[,1],cell_indices[,2])
    return(convergence_measure)
}

PartitionBins <- function(sce, increment_size = 500, seed = 123){
    bins <- c()
    while(dim(sce)[2]>=increment_size){
        if(dim(sce)[2]<increment_size*2){
            bins <- c(bins, sce)
            sce <- sce[,1]
        }
        else{
            set.seed(seed)
            train_ind <- sample(seq_len(ncol(sce)), size = increment_size)
            bins <- c(bins, sce[,train_ind])
            sce <- sce[,-train_ind]
            cat(".")
        }
    }
    return(bins)
}

FindOptimalIncrementSize <- function(sce){
    # Find the computationally optimal increment_size for grainy alignments using 5 increment points (100, 500, 1000, 2500)
    res <- c(100,250,500,750,1000,1500,2500,5000)
    runtimes <- c()
    for(i in res){
        runtimes <- c(runtimes, suppressMessages(system.time(slingshot(sce[,0:i])))["elapsed"])
    }
    mat <- as.data.frame(t(rbind(res,runtimes)))
    c <- summary(lm(mat$runtimes ~ poly(mat$res, 2, raw=TRUE)))$coefficients[,"Estimate"]
    predicted_runtimes <- c()
    for(i in seq(from=100, to=5000, by=1)){
        predicted_runtime <- (c[1] + c[2]*i + c[3]*i^2)*(dim(sce)[2]/i)
        predicted_runtimes <- rbind(predicted_runtimes,c(i,predicted_runtime))
        cat(".")
    }
    increment_size <- predicted_runtimes[which.min(predicted_runtimes[,2]),1]
    return(increment_size)
}

```

# Answer
> 1 votes
There is a recent publication that compares methods for trajectory inference. In Figure 2 in terms of scalability PAGA, PAGA Tree and Component 1 get the highest rating. The authors offer a shiny app to help select the best method for your dataset.

---
Tags: r, scrnaseq, statistics
---


# Question
Title: Decontaminating RNA seqs for de novo transcriptome assembly and annotation of novel eukaryotes
I have raw paired-end RNA-seq reads for two novel eukaryotic species. Some background: the reads represent a copepod (arthropod) species each. The mRNA for each read set was obtained by extracting mRNA from multiple individuals. I am not exactly sure what sequencer was used (some Illumina short read machine as far as I can tell; I am not directly involved in this project).

I'd like to assemble these reads into transcriptomes *de novo* as no reference genomes are available as far as I know; my hope is to assemble with `Trinity` and then annotate with `Dammit` using its default databases (and check for completeness w/ `BUSCO`).

However, it looks like something went wrong during the sequencing runs (or during sample preparation) because the sequence GC content plots from FastQC look like this:

Here the curve with the tall peak is one sample, and the more flatter curve is the other. There are two closely related curves per sample because `FastQC` apparently seems to have recognized each read (forward and reverse) as a separate sample.

I first presumed that this is just rRNA contamination (I believe the mRNA was prepared via poly-A selection), and I filtered the reads against all of `sortmerna`'s default databases (which I believe include sequences from `SILVA` and `Rfam` among others). Then I passed the `sortmerna`'d read through `FastP` to filter out adapters, paired read correction, and low complexity filtering (alongside whatever else `FastP` does by default). I had ran `FastQC` before `FastP` to check if the rRNA filtering did anything and it hadn't, so I also asked `FastP` to filter by length (`--length_required 151`) since there were small amounts of reads with lengths shorter than this value which I thought could be the source of the additional peaks in the data.

However, even after `sortmerna` and `FastP` together, nothing has really changed:

(There are only two traces/curves here instead of four because the paired reads were apparently recognized properly by FastQC now.)

My plan was to run `Khmer` to clean the reads a little bit more and then assemble them. I did do a quick assembly with `Trinity` (with the raw reads) and the `BUSCO` scores were pretty bad (only ~ 25% of the arthropod `BUSCO` genes present in the assemblies). I was hoping I'd be able to get `Dammit` to annotate whatever was left after the cleaning, despite the poor BUSCO prognosis.

Now I'm not entirely sure what to do if I can't clean the reads to begin with, since preventing nonsensical contigs would be absolutely important for a decent annotation.

So this leads me to my questions:

* Are these data even salvageable?
* Is there anything else I can do to clean the data?
* Is this even contamination to begin with that I'm seeing here? (I am presuming some kind of contamination because both samples share that peak at ca. 18% GC content, among other shared features.) Would mapping all the reads against each other and discarding the shared sequences help? If yes, is there a tool I can use for this?
* Any other thoughts in this regard?

I searched through the posts here and on Biostars as well as I could but I wasn't really able to find anything helpful in this regard (hence my own post).

I'm very new to all of this, so any inputs, tips, advice, and/or guidance would be most appreciated.

# Answer
according to your Fastqc files, your reads have a moderate quality. But I am sure there is something to do to improve them. Try to run `fastp` with the following options:

```
fastp -i _1.fq -I _2.fq --detect_adapter_for_pe -l 36 -5 -3 -q 20

```

and run `FastQc` again to check to efficacity of the trimming.
> 1 votes

---
Tags: rna-seq, sequence-annotation, assembly, transcriptome, quality-control
---


# Question
Title: How do I convert an Excel file (.xlsx) to an .smi file format?
How do I convert an Excel file(.xlsx) to .smi file format?

The original file I have is an excel file containing over 27,000 ligands. I am going to using the Padel Descriptor software to calculate the descriptors for the ligands however, Padel Descriptor just supports .smi files. Screenshot of the excel file added below for more clarity.

# Answer
> 1 votes
So a smi file is simply a text file formatted `SMILES` tab `name` newline.

Therefore what you need to do is

* save for safety
* make a new sheet
* copy the SMILES column and the name column
* save as tab separated file
* in Finder/Nautilus/File Explorer rename the file extension to `.smi`

---
Tags: proteins, cheminformatics, computational-biochemistry, small-molecules
---


title: classifying samples by TCGA signature

text:
# Question
Title: classifying samples by TCGA signature
I have some RNA-seq samples from multiple glioblastoma tumours that I'm now trying to classify according to a specific gene signature (from Verhaak et al., 2010) using R. The gene signature is reported as a gene list with specific centroids for each of the 4 clusters (https://api.gdc.cancer.gov/data/941f81a1-05d7-4f84-80ec-534b8dc1ebac). I'm wondering how I can use this signature to classify my samples in R? Would it involve some sort of k-nearest neighbours method?

Additionally, the signature was identified using microarray data, but I am classifying RNA-seq data. Is there any sort of adjustment I should make to the signature to account for this?

Thanks in advance!

# Answer
This is not a very simple problem. I saved the table and "test.txt" and visualize it below:

```
library(pheatmap)
x = read.delim("test.txt")

head(x)
  Gene.Symbol            Composite_chr_coords  Proneural      Neural
1      CDKN1B   [36.1:12:12761576-12766570:+] 0.03353670  0.03353670
2       EPB41    [36.1:1:29086215-29319545:+] 0.03882308  0.03882308
3        CLGN  [36.1:4:141529057-141568232:-] 2.12796056  0.41999847
4      PDE10A  [36.1:6:165664525-165995575:-] 1.66012916  0.25349464
5     RALGPS2  [36.1:1:176960923-177152465:+] 0.78971500 -0.01945154
6        TAF5 [36.1:10:105117714-105138812:+] 1.08354870  0.09461930
    Classical Mesenchymal
1  0.03353670 -0.84379400
2 -0.84821990  0.03882308
3  0.41999847  0.41999847
4  0.25349464  0.25349464
5 -0.01945154 -0.01945154
6  0.09461930  0.09461930

```

You can see there are no so called markers genes and if you look at the scale of the data, it was normalized in microarray and it will be very hard for you to get your RNA-seq to fall into that range, and apply the medoids to identify your clusters or type.

Since you did not provide any data, and it is not trivial to simulate data for this, I can suggest the approach in words here:

1. Subset your count table to these exact set of 840 genes.
2. Normalize the data and do a PCA or hierarchical clustering to see if you indeed observe 4 or 3 or 2 subtypes
3. If you indeed see 4 or 3 clusters this would be straight forward, take the centres of means of these clusters, and correlate with each of the sub types here. Chances are the subtype that is most correlated with one cluster, is the subtype you are interested in.
> 1 votes

---
Tags: r, rna-seq, clustering
---


# Question
Title: How to create a list of differentially expressed (DE) genes after normalization with RUVSeq?
I am using edgeR to perform differential expression (DE) analysis on a set of RNA-seq data samples (2 controls; 8 treatments). To correct for batch effects, I am using RUVSeq.

I am able to get a list of DE genes without normalization:

```
x <- as.factor(rep(c("Ctl","Inf"),c(2,8)))
set <- newSeqExpressionSet(as.matrix(counttable),phenoData=data.frame(x,row.names=colnames(counttable)))
design <- model.matrix(~x, data=pData(set))
y <- DGEList(counts=counts(set), group=x)
y <- calcNormFactors(y, method="upperquartile")
y <- estimateGLMCommonDisp(y, design)
y <- estimateGLMTagwiseDisp(y, design)
fit <- glmFit(y, design)
lrt <- glmLRT(fit, coef=2)
top <- topTags(lrt, n=nrow(set))$table
write.table(top, paste(OUT, "DE_genelist.txt", sep=""))

```

Then immediately after creating the "top" object, I use RUVg to normalize:

```
# [...]
top <- topTags(lrt, n=nrow(set))$table
empirical <- rownames(set)[which(!(rownames(set) %in% rownames(top)[1:5000]))]
ruvg <- RUVg(set, empirical, k=1)
write.table(ruvg, paste(OUT, "DE_RUVg_genelist.txt", sep=""))

```

And I get the error:

```
Error in as.data.frame.default(x[[i]], optional = TRUE) :
  cannot coerce class ‘structure("SeqExpressionSet", package = "EDASeq")’ to a data.frame

```

I am not sure how to print the list of normalized results like I can with the unnormalized data. Ideally, I would get a file with the same format as the edgeR output (as a .csv or .txt file):

```
"logFC" "logCPM" "LR" "PValue" "FDR"
"COBLL1" -2.150 4.427061248733 75.0739519350016 4.53408921348828e-18 9.51203608115384e-15
"UBE2D1" -2.178 3.577168782408 74.9346752854903 4.86549160161322e-18 9.51203608115384e-15
"NEK7" -2.404 4.020072739285 72.6539117671717 1.54500340443843e-17 2.71843349010941e-14
"SMC6" -2.300 5.674738981329 61.8130019860261 3.7767230643666e-15 3.4974443325016e-12

```

**How can I get a list of genes as an output *after* normalization with RUVSeq?**

# Answer
> 1 votes
You do the normalization before running your edgeR. The purpose of `RUVg` is to remove "Remove Unwanted Variation Using Control Genes". In your code, you ran edgeR and then normalize the data using RUVg, which is only going to return you the normalized counts.

Using the example dataset in vignette:

```
library(RUVSeq)
library(zebrafishRNASeq)
data(zfGenes)
filter <- apply(zfGenes, 1, function(x) length(x[x>5])>=2)
filtered <- zfGenes[filter,]
genes <- rownames(filtered)[grep("^ENS", rownames(filtered))]
spikes <- rownames(filtered)[grep("^ERCC", rownames(filtered))]

x <- as.factor(rep(c("Ctl", "Trt"), each=3))
set <- newSeqExpressionSet(as.matrix(filtered),
                           phenoData = data.frame(x, row.names=colnames(filtered)))
set <- betweenLaneNormalization(set, which="upper")

set1 <- RUVg(set, spikes, k=1)

```

You can look at it, it's an expression set with counts etc, not results:

```
set1
SeqExpressionSet (storageMode: lockedEnvironment)
assayData: 20865 features, 6 samples
  element names: counts, normalizedCounts, offset
protocolData: none
phenoData
  sampleNames: Ctl1 Ctl3 ... Trt13 (6 total)
  varLabels: x W_1
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation:

```

You run edgeR now on the results of RUVg:

```
design <- model.matrix(~x + W_1, data=pData(set1))
y <- DGEList(counts=counts(set1), group=x)
y <- calcNormFactors(y, method="upperquartile")
y <- estimateGLMCommonDisp(y, design)
y <- estimateGLMTagwiseDisp(y, design)
fit <- glmFit(y, design)
lrt <- glmLRT(fit, coef=2)
topTags(lrt)

```

# Answer
> 0 votes
I have not used this package but from your code it seems that `ruvg` is not a table. Instead, it is an R object, which means that you cannot use `write.table`. I think the results you want is stored in the object. All R objects contain "slots" of data, which can be accessed by `@`. If I were you, I would type `ruvg@` and should be able to see which data slots are contained in the object.

---
Tags: r, rna-seq, differential-expression, normalization, edger
---


# Question
Title: BWA: Detecting Variation between Reference Genome and Short-Read Sequences
I need to identify all loci in the short-read sequence at which the number of microsatellite repeats (i.e. number of copies of "AA," "GTC," etc.) differ from the reference genome, as well as the locations at which the number of repeats is the same in both the reference genome and the short-read sequence.

I used the Burrows-Wheeler Aligner (BWA mem) to map a high-coverage short-read sequence (obtained from the NCBI's short-read sequence archive) to a reference genome. The output is in .sam format. I have also used a separate program to identify the loci in the reference genome at which microsatellites occur.

I would like to identify all the loci in the short-read sequence at which the microsatellite length and loci differ from the reference genome. Does anyone know any tools or packages I could use to read a .sam/.bam file of a short-read sequence mapped to a reference genome and identify specific loci at which the short-read sequence differs from the reference genome? I am using RStudio and have access to my university's supercomputer clusters.

For info on microsatellites, see here: https://en.wikipedia.org/wiki/Microsatellite#:~:text=A%20microsatellite%20is%20a%20tract,locations%20within%20an%20organism's%20genome.

My end goal is to determine the loci at which microsatellite alleles tend to differ within a species, as well as the loci at which microsatellite alleles tend to be identical between all members of the same species.

# Answer
> 3 votes
I haven't done any whole-genome STR analysis from NGS data myself, but are aware of others that have used lobSTR for this. There's also a recent paper \[here\] that compares a few different STR analysis packages (i.e. RepeatSeq, LobSTR, HipSTR, GangSTR). Here's the concluding paragraph:

> In conclusion, all these tools are built to genotype STRs but have different strengths and weaknesses. Based on our analysis there is no clear overall winner. RepeatSeq and HipSTR are the best when considering genotyping error rate even with low coverage. On the other hand, GangSTR has an advantage because it is the only tool among them that can call alleles longer than the read length but shows higher error rate unless looking at only the enclosed class of reads, which in turn would lose the GangSTR's advantage of picking up long genotypes. In addition, GangSTR is the newest tool and so comes with reference files for different reference builds that are periodically updated according to the tool's webpage. The correct choice of a tool and the subsequent filtering depends on the aim of the analysis, and might be influenced by available hardware resources and time limit for running tools.

In other words, without more information about the specific problem you're trying to solve, it's not clear what the best tool to use is (this is a common issue with Bioinformatics problems).

---
Tags: genome, phylogenetics, variation, bwa, reference-genome
---


# Question
Title: Remove Redundant Sequences from FASTA file in Python
I'm attempting to remove redundant sequences from a fasta file (from NCBI). When I execute this code, it returns the number of spots, not the number of sequences. (Number of spots: 408,293, Number of sequences: 40,830)

Please let me know if you have run a code that can remove the duplicate ideas or have any ideas/suggestions.

Thank you SO much.

```
from Bio import SeqIO
import time

start = time.time()

seen = []
records = []

for record in SeqIO.parse("b4r2.fasta", "fasta"):
    if str(record.seq) not in seen:
        seen.append(str(record.seq))
        records.append(record)

#writing to a fasta file
SeqIO.write(records, "no_dupes_b4r2.fasta", "fasta")
end = time.time()

print(f"Run time is {(end- start)/60}")

```

# Answer
Your code (as currently formatted) should generate the correct output, but it is not efficient. Make `seen` a set to improve these O(n) `in` calls to O(1) time complexity- this will speed it up a lot

```
seen = set()
records = []

for record in SeqIO.parse("b4r2.fasta", "fasta"):
    if record.seq not in seen:
        seen.add(record.seq)
        records.append(record)

```
> 5 votes

---
Tags: fasta, python, sequence
---


# Question
Title: Relaxed and sequential Phylip format conversion
I routinely use AlignIO in BioPython for manipulating alignments and in context here moving from fasta format to phylip format. Alot of phylogeny packages accept a relaxed (sequence ID) & sequential (sequence) phylip format. This format is user-friendly, helps double-check for alignment errors and can be the essential input format.

A typical bit of code is,

```
'''
Converts a fasta file to phylip format
'''
from Bio import AlignIO

Class Convert (inheritClass):

    def fa2phy (self):
        fileOut, fileIn = self._inOut() # builds the path name and deletes any pre-exising outfile
        align = AlignIO.read(fileIn, "fasta") # reads the alignment in fasta
        AlignIO.write(align, fileOut, "phylip-relaxed") # writes in phylip

```

The output of "phylip-relaxed" is interleaved format

```
4 120
virusA_detailed_description ACGTGACG .... ATCCCCAA
virusB_detailed_description ACGTGACG .... ATCGGGAA
virusC_detailed_description ACGTGACG .... ATCAAAAA
virusD_detailed_description ACGTGACG .... ATTTTTAA

ACGTGACG .... ATCCCCAA
ACGTGACG .... ATCGGGAA
ACGTGACG .... ATCAAAAA
ACGTGACG .... ATTTTTAA

```

Output sought is relaxed-sequential

```
4 120
virusA_detailed_description ACGTGACG .... ATCCCCAAACGTGACG .... ATCCCCAA
virusB_detailed_description ACGTGACG .... ATCGGGAAACGTGACG .... ATCGGGAA
virusC_detailed_description ACGTGACG .... ATCAAAAAACGTGACG .... ATCAAAAA
virusD_detailed_description ACGTGACG .... ATTTTTAAACGTGACG .... ATTTTTAA

```

There is an option to allow a sequential format (above). HOWEVER, a `"phylip-sequential"` option output is as follows (it truncates the sequence ID to 10 characters,

```
4 120
virusA_det ACGTGACG .... ATCCCCAAACGTGACG .... ATCCCCAA
virusB_det ACGTGACG .... ATCGGGAAACGTGACG .... ATCGGGAA
virusC_det ACGTGACG .... ATCAAAAAACGTGACG .... ATCAAAAA
virusD_det ACGTGACG .... ATTTTTAAACGTGACG .... ATTTTTAA

```

No matter how I look through AlignIO options the above is the nearest format. Could anyone point to a IO style solution? I could write it, but that defeats the object (i.e. purpose not an instance of a method).

---

*Note* the `"...."` in the outputs above are not a typical "sequence identity" symbols, it just means that total length of a line in an interleaved format is 60 nucleotides (it might be 80).

# Answer
I have nutted out a solution for you.

If this is the content of input.fa:

```
>virusA_detailed_description
ACGTGACGATCCCCAAACGTGACGATCCCCAA
>virusB_detailed_description
ACGTGACGATCGGGAAACGTGACGATCGGGAA
>virusC_detailed_description
ACGTGACGATCAAAAAACGTGACGATCAAAAA
>virusD_detailed_description
ACGTGACGATTTTTAAACGTGACGATTTTTAA

```

And this is the script:

```
#!/usr/bin/env python3

'''
Converts a fasta file to phylip format
'''

from Bio import AlignIO
from Bio.AlignIO.PhylipIO import PhylipWriter
import sys

infile = sys.argv[1]
outfile = sys.argv[2]
idwidth = int(sys.argv[3])
with open(infile, 'r') as input_handle:
    align = AlignIO.read(input_handle, "fasta")
    with open(outfile, 'w') as output_handle:
        PhylipWriter(output_handle).write_alignment(align, id_width=idwidth)

```

Then, save this script as `phylpwrtr.py` and do:

```
chmod +x phylpwrtr.py
./phylpwrtr.py input.fa output_longids.phy 30 # 30 is the width of the header, adjust the value to adjust header length

```

The contents of `output_longids.phy` looks like:

```
 4 32
virusA_detailed_description    ACGTGACGAT CCCCAAACGT GACGATCCCC AA
virusB_detailed_description    ACGTGACGAT CGGGAAACGT GACGATCGGG AA
virusC_detailed_description    ACGTGACGAT CAAAAAACGT GACGATCAAA AA
virusD_detailed_description    ACGTGACGAT TTTTAAACGT GACGATTTTT AA

```
> 1 votes

---
Tags: phylogenetics, python, format-conversion
---


# Question
Title: Is there a way to customize what observations are plotted in ComplexHeatmap?
So right now my observations in the dataset I've processed is stored as sample identifier columns and gene rows: an example of this would be for sample 1 and gene TRIM21 the observation is Missense Mutation;Frame Shift Deletion. I'm having issues separating the two alterations and instead of being able to show the gene in the sample having 2 different alterations, I end up having to generate a new colour classification (being that the gene in the specific sample exhibits specifically both a missense mutation and a frame shift deletion instead of being able to represent the observation with just colours and missense and frame shift deletion colours alone). Here is an a screenshot of what my matrix looks like:

As an example here is some test code generating the test plot:

```
library(ComplexHeatmap)

test_labels <- c("Missense_Mutation", "Splice_Site", "Frame_Shift_Del",
                 "Splice_Region","Missense_Mutation;Splice_Site", "Splice_Site;Frame_Shift_Del",
                 "Splice_Site;Missense_Mutation", "Frame_Shift_Del;Splice_Site")

plot <- Heatmap(test, row_names_side = "left",
        row_names_gp = gpar(cex = 0.7),
        column_names_gp = gpar(cex = 0.7),
        heatmap_legend_param = list(title = "Alterations",
                                    at = test_labels, labels = test_labels), show_heatmap_legend = TRUE)

plot

```

I've attached a screenshot of my actual plot and as you will be able to see, this ends up creating multiple different classifications:

# Answer
> 1 votes
What is fed into `ComplexHeatmap` is what it outputs! So if one of your genes has a `Missense_Mutation` as well as a `Frame_Shift_Del` and you would like to plot its value, of course you will have a category and color that is different than those of `Missense_Mutation` or `Frame_Shift_Del`, that is only natural.

What you can do is to create additional variables like `your_gene_Missense_Mutation` or `your_gene_Frame_Shift_Del` each with 0 or 1 to indicate if a particular mutation type is present or not for a particular gene (you can think of it like one-hot-encoding). With a lot of different levels, the resulting heatmap could be difficult to interpret though.

Another simpler alternative would be a using a "universal label" if your gene has multiple mutations. Basically if the cell contains a `;`, its value will be changed to something like `mixed_mutation_status` and instead of `multiple "mutation_type_1; mutation_type_2; ...` like categories, you will just have one `mixed_mutation_status` on top of `Missense_Mutation`, `Frame_Shift_Del`, ...

---
Tags: bioconductor, mutations, heatmap, complexheatmap
---


# Question
Title: Is there an "official" standard usage of Red and Green on differential gene expression heatmaps?
I have recently found myself making multiple heatmaps for visualization of differential gene expression results between replicates in two experimental conditions.

The default settings in the plotting software I use color:

* Upregulated genes --\> red
* Downregulated genes --\> green
* Unchanged genes --\> black

My intuition *screams* at me that this color scheme is reversed, and that upregulated genes should be green and that downregulated genes should be red. I have asked several labmates about their intuition regarding this color scheme, and have found about a 50/50 split in those that intuitively agree with each different scheme.

In my searching on this topic, I have found that several programs plot heatmaps with the same default colors as the program I'm using (red=upregulated, green=downregulated). Additionally, I have been told in conversation that this is essentially an accepted standard.

If this is indeed the convention, I will begin trying to rewire my intuition, but I want to confirm that it is before I do so. However I have not yet been able to find any "official" guidelines from agencies, publishers, authorities in the field, etc., specifying the use of red/green in heatmaps.

I do know that there are many suggestions for alternative heatmap color schemes, but I'm interested in finding an official answer for the red/green scheme specifically.

**Edit:** Several comments have stated that there is no such standard as described above. Therefore I am modifying the following question somewhat to:

**Are there "official" guidelines from agencies, publishers, authorities in the field, etc., specifying the use of red/green in heatmaps? Or conversely are there any “official” guidelines stating that there is no such specification and/or that any such selection of colors is acceptable?**

# Answer
> 2 votes
As others already pointed out, there is no official heatmap color scheme. I am not sure there is an official color scheme for any other types of plots either.

However, heatmaps seem to be far more standardized. Traditionally, heatmaps have been red-black-green as you originally stated. Most heatmaps from the microarray era and early RNA-seq are like that. You can have a heatmap for any type of data, but the majority is for gene expression. A lot of tools and tutorials were written back then, so they have that color scheme. From about 2008 to 2012, most people switched to the red-white-blue (or red-yellow-blue) color scheme. Switching green to blue avoids the color-blindness issue. Switching black to white is better for the ink if you print your papers and generally looks better on white background.

I originally found red as high to be confusing because I was thinking in financial terms (if the market is red, it's down). Now I think of it as a water faucet (red is hot, blue is cold).

Edit: In the last year or so, the viridis palette started to become widely used. It's far from the standard, but could become more popular in the future.

# Answer
> 1 votes
Although there is no standard as DeepTools' Devon Ryan said, there are still some habits. Usually red means upregulated, as you can see in the default color scheme in MATLAB's clustergram function. But you also want to make sure color-blinded people can see the difference. So another color scheme might be better, such as red=high while blue=low. Anyway, if you can easily flip the color, using red to represent high is expected.

# Answer
> 0 votes
In the absence of any yet-identified official recommendations, I will update this answer with resources as I find them:

**Example usage in tutorials:**

**Program Defaults:**

**Guides:**

**Tools (by year):**

---
Tags: differential-expression, visualization, heatmap
---


title: Blast results filters

text:
# Question
Title: Blast results filters
I performed a blastn search of NGS data against ssRNA database download from Internet, with a expected value 10-4. The size of NGS data reads is of 125 bp.

I have analyzed the blast results of the libraries. These libraries are from different habitats. The most difficult habitat, oakwood, produced a blast results with a maximum length of 51 bp. These results have a qcov of 100 (in most of the cases).

I would like to know if these low length sizes comes from degradation?

Should do I filter these samples with a length = 125 bp with the aim to obtain results with better quality or with a qcov of 100 is enought to get reliable results?

Thanks you in advance.

# Answer
> 1 votes
I had a think about this one and am stumped. I thought there could be a connection with acute oak disease (AOD) or else chronic oak disease (COD), which are emergent infections threatening the population. However, they appear to be fungal and bacterial infections rather than viral.

Single stranded RNA plant viruses are very common, more so than in mammals, but to achieve the results you are obtaining the trees would need to be really sick if the predominent RNA was viral. Short viral RNAs are produced regularly by viruses as defective RNA and its a major area of research, considering a mechanism of defective interference,

In context it would be weird if this RNA swamped the tree/plant, you need an actively replicating virus for this stuff to appear. So I don't think its either viral or defective viral RNA.

The antiviral response in mammals can involve RNase and this could explain your degradation. Essentially, mammals (probably vertebrates) kick-back hard when virus appears. However, I couldn't really find evidence of this in plants.

I don't think the results you obtained have a clear biological explanation in pathogenesis. The best explanation is it was wet lab-technique, e.g. library construction resulted in overdigestion/fragmentation. The only thing I can suggest is template assembly using a normally distributed plant/oak RNA library.

I was probably overthinking in forging a link with an emergent (tree) pathogen.

---

!

---
Tags: blast, ngs, blastn
---


# Question
Title: What does aov(glm model) actually do in R?
What is happening when you use `aov()` on a glm model in R?

Normal glm model `sumamry(m1)`:

```
Call:
glm(formula = value ~ origin + variable, family = poisson(),
    data = comb_data)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-0.8227  -0.5547  -0.4385  -0.3038   2.5258

Coefficients: (1 not defined because of singularities)
                                                                  Estimate Std. Error z value Pr(>|z|)
(Intercept)                                                     -2.853e+00  5.774e-01  -4.941 7.78e-07 ***
originisolate                                                    4.700e-01  7.071e-01   0.665  0.50625
variableRB17                                                     2.877e-01  5.401e-01   0.533  0.59425
variableM56                                                      1.299e+00  4.606e-01   2.821  0.00479 **
variableRB68                                                    -6.269e-14  5.774e-01   0.000  1.00000
variableM39b4                                                   -1.010e-13  5.774e-01   0.000  1.00000
variableRB108                                                    1.542e-01  5.563e-01   0.277  0.78172
variableM32                                                      4.055e-01  5.270e-01   0.769  0.44171
variableRB33                                                    -6.931e-01  7.071e-01  -0.980  0.32696
variableM41                                                      6.061e-01  5.075e-01   1.194  0.23236
variableRB1                                                      6.931e-01  5.000e-01   1.386  0.16566
variableM39w2                                                    1.542e-01  5.563e-01   0.277  0.78172
variableRB20                                                    -1.792e+00  1.080e+00  -1.659  0.09714 .
variableRB99                                                     1.542e-01  5.563e-01   0.277  0.78172
variableRB110                                                    6.931e-01  5.000e-01   1.386  0.16566
variableRB5                                                      1.542e-01  5.563e-01   0.277  0.78172
variableRB29                                                    -1.823e-01  6.055e-01  -0.301  0.76334
variableRB56                                                    -7.411e-14  5.774e-01   0.000  1.00000
variableRB24                                                    -1.099e+00  8.165e-01  -1.346  0.17846
variableStreptomyces_alboflavus_strain_MDJK44_RB5                1.204e+00  6.583e-01   1.829  0.06740 .
variableStreptomyces_violaceoruber_RB1                           1.299e+00  6.513e-01   1.995  0.04607 *
variableStreptomyces_T44_RB17                                    1.204e+00  6.583e-01   1.829  0.06740 .
variableStreptomyces_violaceoruber_RB1_RB110                     1.299e+00  6.513e-01   1.995  0.04607 *
variableActinomadura_verrucosospora_RB99                         5.108e-01  7.303e-01   0.699  0.48425
variableStreptomyces_autolyticus_strain_CGMCC0516_M56            1.735e+00  6.262e-01   2.770  0.00561 **
variableStreptomyces_atratus_RB108                               9.808e-01  6.770e-01   1.449  0.14740
variableStreptomyces_sp_T44_RB17                                 1.204e+00  6.583e-01   1.829  0.06740 .
variableStreptomyces_sp_QMT.28_M41                               8.473e-01  6.901e-01   1.228  0.21950
variableAmycolatopsis_sp_AA4_M39                                 8.473e-01  6.901e-01   1.228  0.21950
variableLuteimicrobium_xylanilyticum_RB22                       -1.099e+00  1.155e+00  -0.951  0.34138
variableMycolicibacterium_parafortuitum_RB24                     1.239e-13  8.165e-01   0.000  1.00000
variableNocardia_nova_SH22a_RB20_RB56                            2.877e-01  7.638e-01   0.377  0.70642
variableActinomadura_atramentaria_RB29                           5.108e-01  7.303e-01   0.699  0.48425
variableMycolicibacterium_madagascariense_strain_JCM_13574_RB33         NA         NA      NA       NA
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 1024.95  on 1949  degrees of freedom
Residual deviance:  943.36  on 1917  degrees of freedom
AIC: 1507.4

Number of Fisher Scoring iterations: 6

```

A `summary(aov(m1))`:

```
              Df Sum Sq Mean Sq F value   Pr(>F)
origin         1   0.33  0.3285   3.050   0.0809 .
variable      31  10.41  0.3357   3.117 2.06e-08 ***
Residuals   1917 206.47  0.1077

```

# Answer
It's not appropriate to run `aov` on a GLM model. `aov` assumes normally distributed data and it runs on a GLM as if it was a least squares linear model. It uses the working responses and working weights from the GLM fit and performs an ordinary anova as if the responses were normally distributed.

On the other hand, `anova` is a generic function with a method for GLMs and it performs a sequential analysis of deviance.

Purely statistical questions like this are best sent to CrossValidated, although they don't like questions that are posed entirely in terms of an R function. Questions purely about the functionality of functions in the core distribution of R are best sent to R-help, which is a mailing list rather than web forum (see https://www.r-project.org/mail.html and https://stat.ethz.ch/mailman/listinfo/r-help).
> 4 votes

# Answer
It's equivalent to doing `anova(m1, test="F")`

You probably want to do some sort of deviance test for a Poisson glm

```
anova(m1, test="Chisq")

```
> 0 votes

---
Tags: r, statistics
---


# Question
Title: how to have a blast xml file in a readable and understandable way
I got a blast xml file after a DNA sequence search. The XML file is not very readable. What is the way to make it more understandable? This blast contains about twenty matches. Here is the xml file

# Answer
The most human-readable blast output formats are 0-4, e.g. `-outfmt 0` \- which if I am not mistaken is the format used by web blast. However, these formats are a pain to automatically parse.

As Cowboy\_Patrick pointed out, XML is more of a machine-readable format. The tabular output format ( `-outfmt 6`) is very commonly used, because it is easily parsed by unix command-line tools, and its structure also allows a somewhat easy human reading and understanding.
> 2 votes

# Answer
Is this what you are looking for: https://codebeautify.org/xmlviewer ? XML is a format designed to be readable firstly by machines, not humans.

With xslt you can convert XML to various more readable formats see this blog for something that fits your XML: http://plindenbaum.blogspot.com/2008/05/ncbi-blast-xslt-xhtml-svg.html?m=1
> 1 votes

---
Tags: blast
---


# Question
Title: Comparing gene abundances between metagenomes
My workflow until now:

Find fragments of a marker gene in unassembled metagenomes \> download and assemble metagenomes \> recover the gene neighborhood / gene set of interest

Right now I have a rough estimate of how abundant these genes are by the 'depth' noted on the assembled contigs (I use MEGAHIT for assembly). I was wondering if there's a more thorough/proper way to do this. I would like to compare the abundance of specific genes between a) samples in the same study, and b) different studies. I imagine that the size of the individual metagenomes should be considered in both cases, but point b) might add additional difficulties such as different sequencing techniques. I would appreciate your insights.

# Answer
> 2 votes
I would avoid using assemblies to answer this question, as there's no guarantee that you will be able to assemble your genes of interest; you can however estimate their abundance even if they are relatively rare.

How I understand your question as being one of estimating the abundance of either some specific genes (e.g. butyrate metabolism genes) or *all* genes in a microbial community across multiple samples for comparative purposes. In other words, not 16S or marker gene analysis for the purposes of estimating organismal abundance, which is a rather different problem (though in that case I would still not use an assembly).

A more standard workflow is:

1. align metagenomic reads against some existing database of genes annotated appropriately.
2. estimate the number of reads aligning against some gene or ortholog in some way (using for example KEGG Orthology or similar).
3. use counts from (2) as input to some statistical procedure, possibly summarizing across functional categories.

Some examples of how this has been done are here, here, here. I am sure that there are more recent/relevant references but I haven't been following the field closely in the last few years.

---
Tags: gene, assembly, metagenome
---


title: Mapping UniProt id to Entrez id

text:
# Question
Title: Mapping UniProt id to Entrez id
How can we get mapping from UniProt id to Entrez ID? I have UniProt id, but want to convert those ids to Entrez or Ensembl id.

# Answer
The Retrieve/ID mapping tool from Uniprot allows you to convert from Uniprot to entrez ID and vice-versa. It is available here. Just select the appropriate values in the `2.Selection options` panel.
> 2 votes

---
Tags: r, identifiers
---


# Question
Title: Bash script for automatic overlap merging of respective paired-end reads from samples
I have around 300 paired-end files from 150 samples (1 forward and 1 reverse read for each sample). I want to merge (by sequence overlap) respective forward and reverse reads for each of the samples. The samples are named e.g. **SRR2155174\_1.fastq SRR2155174\_2.fastq SRR2155319\_1.fastq SRR2155319\_2.fastq.** Can anyone please help me write a bash script to do this?

# Answer
> 1 votes
Your original question mentioned concatenating things, which is not what you want (that would result in all of the reads from the \_2 files being after all of those in \_1). What you want is to merge overlapping paired end reads. There are a variety of tools available for that, from flash to bbmerge. The general structure then is:

```
for R1 in *_1.fastq ; do
    sampleName=${R1%%_1.fastq}
    some command with ${sampleName}_1.fastq and ${sampleName}_2.fastq producing ${sampleName}.fastq
done

```

You can figure out the `some command ...` portion yourself.

---
Tags: fastq
---


# Question
Title: Python: How to write duplicate sequences removed from fasta file to new file
I currently am using this code to remove duplicate sequences from the fasta file. However, I would also like to write a new file with only the removed duplicates as well as a count for how many times they appear. Please let me know if you have any suggestions, I am hoping to keep it written in one py file, but I am not opposed to running a second file to do so. I attended to write the appended strings to a new file but it looked slightly off...

from Bio import SeqIO
import time

start = time.time()

seen = []
records = []

 for record in SeqIO.parse("b4r2.fasta", "fasta"):
   if str(record.seq) not in seen:
       seen.append(str(record.seq))
   records.append(record)

#writing to a fasta file
SeqIO.write(records, "no_dupes_b4r2.fasta", "fasta")
end = time.time()

print(f"Run time is {(end- start)/60}")
```
</code></pre>

# Answer
> 3 votes
**Edit: user @RamRS has a more robust answer for actually using this. Please use that. Leaving answer up as OP appears to be interested in implementing as a learning activity.**

**Please let us know why this is not doing what you intend.**

I haven't checked it, but something like this would work:

```
from Bio import SeqIO
import time
from collections import defaultdict

start = time.time()

seq_counts = defaultdict(0)
records = []

 for record in SeqIO.parse("b4r2.fasta", "fasta"):
   seq = str(record.seq)
   if seq in seq_counts:  # python is smart about dict keys
       records.append(record)

   seq_counts[seq] += 1

#writing to a fasta file
SeqIO.write(records, "no_dupes_b4r2.fasta", "fasta")
end = time.time()

# print the counts
for seq in seq_counts:
    print(seq_counts[seq], seq)

print(f"Run time is {(end- start)/60}")

```

It will be a bit slow for large fasta files and I am sure that some people will be very upset that I'm using a defaultdict for this.

The code that you have would probably not work as you would want the `records.append()` call to be only run in the if statement testing for non-membership.

There is absolutely no reason to have to write more than one script for this.

# Answer
> 4 votes
Please do not reinvent the wheel - `seqkit rmdup` does what you want.

From the website:

---

# rmdup

## Usage

```
remove duplicated sequences by id/name/sequence

Usage:
  seqkit rmdup [flags]

Flags:
  -n, --by-name                by full name instead of just id
  -s, --by-seq                 by seq
  -D, --dup-num-file string    file to save number and list of duplicated seqs
  -d, --dup-seqs-file string   file to save duplicated seqs
  -h, --help                   help for rmdup
  -i, --ignore-case            ignore case

```

How does a `-D` file look like, you ask:

```
$ zcat hairpin.fa.gz \
    | seqkit rmdup -s -i -o clean.fa.gz -d duplicated.fa.gz -D duplicated.detail.txt

$ cat duplicated.detail.txt   # here is not the entire list
3   hsa-mir-424, mml-mir-424, ppy-mir-424
3   hsa-mir-342, mml-mir-342, ppy-mir-342
2   ngi-mir-932, nlo-mir-932
2   ssc-mir-9784-1, ssc-mir-9784-2

```

---
Tags: fasta, python
---


# Question
Title: RNA seq .counts.txt to bigwig conversion
I've downloaded neuronal RNA seq data from GEO. The files are in .counts.txt format. Could I convert them to bigwig? If so, how?

# Answer
> 1 votes
No, count-level data doesn't have the base-level resolution that is required in order to be usefully visualised on a genome browser. You need a BAM file (or raw FASTQ files that can be used to generate a BAM file) in order to do that.

"No" isn't completely correct: it'd be *possible* to take gene counts and information about gene annotation boundaries to shoehorn them into a BedGraph / BigWig file, but whatever comes out the other end of that isn't going to be a good representation of reality.

---
Tags: rna-seq, rna, bigwig
---


title: Selecting part of an extracted ligand

text:
# Question
Title: Selecting part of an extracted ligand
I have the PDB 6LU7.

* I extracted the ligand using extract lig, org and then the protein using extract prot, poly. The ligand has three parts: 02J, PJE, and 010.
* I want to select the 02J and 010 IN A SCRIPT and delete them.

I know how to do it by simply clicking, but want to incluse this in a script.

# Answer
I know that structure way too well —SARS-COV-2 MPro.

In PyMOL command line

```
 remove resn 02J or 010

```

or

```
 cmd.remove('resn 02J or 010')

```

In python

```
 import pymol2
 with pymol2.PyMOL() as pymol:
      pymol.cmd.fetch('6LU7')
      pymol.cmd.remove('resn 02J or 010')
      pymol.cmd.save('6LU7.mod.pdb')
      #or pdbblock = pymol.cmd.get_pdbstr()

```

From the sounds of it you'd like to make the peptidomimetic a peptide. If you wanted to do the horrid thing of making PJE, the capped side chain glutamine ketamide, into glutamine, the command to rename atoms is a multitude of `cmd.alter('resn PJE and name <atomname>', 'name="<newname>")` followed by a `cmd.alter('resn PJE', 'resn="GLN"')` and lastly `cmd.sort()` to commit the changes. Not pretty, but is all in PyMOL.
> 1 votes

---
Tags: covid-19, sars-cov-2, pymol
---


# Question
Title: Programmatically retrieve Accession List from NCBI
I would like to programmatically retrieve the file "Accession List" from an NCBI page.

I'm running firefox, so I'm clicking the file and clicking `Ctrl` \+ `Shift` \+ `E` to try and figure out what the file endpoint is. But I can't seem to figure it out (I haven't done this much, I just remember doing it once for a previous project).

When I do this, I see two links popping up on my screen, as shown below:

I see that the first one has a field (I'm talking about the panel on the right) called "acc", which varies from one data dump to the next, so I'm presuming this is the link I need to use to retrieve the file (though I'm not sure), and just change the "acc" field depending on the data I want. Assuming this is right, I'm not sure how to go from here. I've included the text on the right below for clarity.

```
GET
    https://www.google-analytics.com/collect?v=1&_v=j83&a=1402353855&t=event&_s=4&dl=https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA491191&o=acc_s%3Aa&ul=en-us&de=UTF-8&dt=Run Selector :: NCBI&sd=24-bit&sr=1280x720&vp=1232x601&je=0&ec=Link Click&ea=Metadata  &el=NONE&_u=SACAAEAB~&jid=&gjid=&cid=1350895533.1590937485&tid=UA-68648221-3&_gid=391407234.1594002343&cd1=run_selector&cd2=sra&cd3=&cd5=NO&cd6=run_selector:sra:&cd8=&cd9=&cd10=&cd11=&cd12=&cd13=&cd14=&cd15=&cd16=&cg1=run_selector&cg2=sra&cg3=&cg4=run_selector:sra:&z=508671370

```

# Answer
Which accessions are you trying to get to? In the page you refer to, there is an option to download the list of accessions - one SRR run accession per line - by clicking the 'Accessions' button or the entire table by clicking the 'Metadata' button as highlighted in the screenshot below:

You can get the same data from the command line using the Entrez Direct. The command for that would be something like:

```
esearch -db sra -query 'PRJNA491191[bioproject]' | efetch -format runinfo

```

to download the table in CSV format.
> 2 votes

---
Tags: phylogenetics
---


# Question
Title: GSEA enrichr with 10x genomics differential_expression ranks
I am attempting to use GSEA enrichr with 10x genomics differential\_expression rankings. Reading what people seem to be doing with GSEA, there seems to be a pre-/post-singlecell gene expression versioning of the tools associated, so I want to make sure what I am doing makes sense.

Explanation:

For each gene in columns 1 and 2, there are 3 columns in each of the clusters that the 10X Genomics `cellranger count` software produces, e.g.:

```
Cluster 3 Mean Counts   Cluster 3 Log2 fold change   Cluster 3 Adjusted p value
0.2392551160479023      -6.591931364165717           0.029177998608334674
0.1435530696287414      -5.608385347948343           0.0003833420249831818
1.107409394278862       -5.441829814465051           3.8059567638748296e-09
0.1948220230675776      -4.9155704967404             0.008739337761029243
0.20849374398460058     -4.65702722396748            1.0980978990257375e-11
0.36571853453036496     -4.6498340596466585          1.1896434946397147e-17
0.1811503021505546      -4.445327449243366           1.3416233207430882e-12
0.4409129995739914      -4.433552335619653           1.2400396349891802e-16
1.7226368355448967      -4.240367394494225           2.2072280062348652e-15
0.10937376733618391     -3.902215199856403           2.0315030904651526e-12
3.2948847410025404      -3.7580328996321235          0.028858857004902847
0.3588826740718535      -3.738237134795285           8.439629283346646e-13
[...]

```

In trying to give a tissue / cell type identity to each of these clusters, I've used the following method to feed a subset of genes to EnrichR:

1- Sort the `Log Fold change` column numerically so we get low negative values. 2- Filter out for the `Adjusted p-value` to have a value smaller than 0.05. 3- Filter out for the `Means Counts` value to be bigger than 0.01. 4- Take the top n genes from the list, and feed it to GSEA EnrichR to compare against Mouse\_Gene\_Atlas or ARCHS4\_Tissues.

Given that the fold change of the clusters is relative to each other, does it make more sense to rank the top n genes numerically or in inverse numerical value? E.g. selecting the first 100 genes in each cluster that have the most negative `Log2 fold change` values or the most positive `Log2 fold change` values or both (the top 100 genes with highest absolute `Log2 fold change`)?

Looking at previous answers, I get that one shouldn't rank by p-value: e.g. see @gringer reply https://bioinformatics.stackexchange.com/a/13465/180

# Answer
You should be using the most positive fold changes, not the most negative. The positive ones are the ones that are most highly expressed in your cluster of interest. You probably want to use less than 100 genes. Depending on the experiment, you may not have that many good markers for a particular cluster.

The other question you mentioned refers to DESeq2 specifically, so it is not applicable here.
> 1 votes

---
Tags: single-cell, differential-expression, gsea
---