# Question
Title: Converting a data frame to multiple lists
I have data frame each row is a pathway and it involved genes like

like

```
> head(h)
    Release.of.Cancer.Cell.Antigens.1  ARID1A    ATM  BIRC3   BLM BNIP3   BRCA1   BRCA2 BRD3  BRD4 BRIP1 CCNA1 CCND1
V1        Cancer Antigen Presentation    ARG2   ATF3    B2M BATF3   BLK   CCL22 CCL3/L1 CCL4  CCL5  CCR5  CD14  CD19
V2      T-cell Priming and Activation ADORA2A  BATF3   BTLA  CCL2 CCL22 CCL3/L1    CCL4 CCL5  CCR2  CCR5  CD1C   CD2
V3 Immune Cell Localization to Tumors     A2M ADGRE1 ANGPT2  APOE  ATF3     AXL   BATF3  BLK CCL13 CCL14 CCL18 CCL19

```

I want to convert that like

```
gene.sets <- list(Release.of.Cancer.Cell.Antigens = c("ARID1A", "ATM",  "BIRC3",    "BLM"),
            Cancer_Antigen_Presentation = c("ARG2","ATF3","B2M"))

```

I have not done anything so far as I am not certain what exactly I should google for.

By code the user kindly provided I get this

```
gene.sets = lapply(1:nrow(h), function(x) unname(unlist(h[x,-1])))
names(gene.sets) = h[,1]

[91] "IGF2R"     "IHH"       "IL11"      "IL11RA"    "IL12RB2"   "IL15"      "IL17A"     "IL1A"      "IL1B"
[100] "IL22RA1"   "IL2RA"     "IL2RB"     "IL6"       "IL7R"      "IRF1"      "IRF2"      "IRF3"      "IRF4"
[109] "IRF7"      "IRF9"      "ISG15"     "ITGAE"     "JAK1"      "JAK2"      "JAK3"      "KIT"       "LAG3"
[118] "LAIR1"     "LIF"       "LTB"       "MX1"       "NCAM1"     "NKG7"      "OAS1"      "OAS2"      "OAS3"
[127] "PDCD1"     "PDCD1LG2"  "PF4"       "PIK3CA"    "PIK3CG"    "PIK3R5"    "PRF1"      "PRLR"      "PROM1"
[136] "PSMB10"    "PSMB5"     "PSMB9"     "PTEN"      "REN"       "RORC"      "RSAD2"     "SELL"      "SELP"
[145] "SIRPA"     "SOX11"     "SOX2"      "SPRY4"     "STAT1"     "STAT2"     "STAT4"     "TBX21"     "TIGIT"
[154] "TNF"       "TNFAIP3"   "TNFRSF10B" "TNFRSF10C" "TNFRSF11A" "TNFRSF11B" "TNFRSF14"  "TNFRSF17"  "TNFRSF18"
[163] "TNFRSF1A"  "TNFRSF1B"  "TNFRSF4"   "TNFRSF8"   "TNFRSF9"   "TNFSF10"   "TNFSF12"   "TNFSF13"   "TNFSF13B"
[172] "TNFSF18"   "TNFSF4"    "TNFSF8"    "TPSAB1/B2" "TRAF1"     "TSLP"      "VEGFA"     "WNT7B"     ""
[208] ""          ""          ""          ""          ""          ""          ""          ""          ""
[217] ""          ""          ""          ""          ""          ""          ""          ""          ""
[226] ""          ""          ""          ""          ""          ""          ""          ""          ""
[235] ""          ""          ""          ""          ""          ""          ""          ""          ""
[244] ""          ""          ""          ""          ""          ""          ""          ""          ""
[253] ""          ""          ""          ""          ""          ""          ""          ""          ""

x=lapply(gene.sets, function(z){ z[!is.na(z) & z != ""]})

```

# Answer
> 2 votes
I don't know whether this will actually work since you didn't provide a reproducible example...

```
gene.sets = lapply(1:nrow(h), function(x) unname(unlist(h[x,-1])))
names(gene.sets) = h[,1]

```

---
Tags: r
---


# Question
Title: Extract compounds/proteins from plant genomic data
I want a process by which I can identify genes within plant genomes; i.e. a process that works for all plants, not just for a single species.

I'm wondering, if, by bioinformatic means only, it is possible to "extract" genes and, after that proteins and natural compounds from plant genome data that is available on NCBI (please see https://www.ncbi.nlm.nih.gov/datasets/genome/GCA\_029618835.1/ as an example)? Example: just by processing the data from the link above we can get to Flavoxanthin or any other substance.

What other useful information from a chemical compound perspective can be found from digging into the genome itself?

Any hints will be greatly appreciated.

Thank you!

# Answer
I believe I understand the project, to obtain a biochemical/metabolic profile of a large number of poorly annotated plant genomes. There isn't a perfect solution, like Prokka for bacteria. It is clearly an area subject to recent development and interest.

I would recommend using Gene Ontology Meta Annotator for Plants (GOMAP) which is described in "Plant Methods" journal.

There are limitations:

* KEGG is not used (below);
* its a singularity container - and thats its only distribution.

**Learning Singularlity** You will need to learn running the container through Singularity and that will be limiting if you are not familiar with containers. However, this is the best of the available "out of the box" solutions I can see in plants. A Singularity HPC is a good call here, tricky if you don't have access to one. I know Docker, Singularity is similar but different. The information on Singularity is here. I would start reading from "Executing Commands"

```
singularity pull docker://godlovedc/lolcow # get the container
singularity exec lolcow_latest.sif cowsay moo # run it

```

(the cow is a famous teaching tool for containers)

**KEGG** One approach that can be used is to perform COG (NCBI/Genbank code) and KEGG no. identification (e.g. eggNOG). Tools such as eggNOG are setup for metagenomics data however. The COG/KEGG no. output is only the starting point, interrogating each of these to assign gene function and metabolic pathways will require a KEGG database and custom code. Unless someone has published, e.g. KEGG database in SQL.

In my experience the Kyoto Encyclopaedia of Genes and Genomes (KEGG) is the best singular functional database available, thats why I see GOMAPs exclusion of this as a limitation. COG doesn't provide a particularly useful information functional description IMO. eggNOG, or any other method to identify COGs and KEGG numbers, interrogated against the KEGG database works for plants very well. I know this because I'm not focused on plants but often recover plant defence mechanisms genes from KEGG. Basically, KEGG identifies \>\> half of all genes with known COG nomenclature and notably higher than other functional databases.Thus I see GOMAPs exclusion of KEGG as a limitation. However, to reiterate to go identify COG/KEGG -\> KEGG database interrogation and analysis is a lot of custom code.

**Databases** GOMAP implements *Arabidopsis* annotation, which may be cool and provides a starting point. It's also running UniProt, UniProt-plants, interpro, GO databases and those **+** *Arabidopsis* annotated DB could be really cool. It definitely sounds promising.

**Processor time** Any method of gene identification takes AGES and I have no doubt that GOMAP will be REALLY slow, looking at the pipeline it's using - there are a lot of databases its using. This is GOMAPs justification for using Singularity, it going to take ages to run therefore containerise it in a format often used by modern HPCs. To be honest it's a bit lazy, especially as it was published in 2021 (should be uploaded e.g. to conda). The authors talk about Python and R ... hmmmm, it looks to be a MySQL backend, thats where its speed will come from. The container will be large BTW hopefully is \<1GB, but it will be \>200MB.

**Summary** Its worth learning Singularity just keep to the basic commands to run GOMAP. You will need access to a large machine and the output will be okay but they're better *general* databases. However, I don't know how useful the *Arabidopsis* database it - thats about the biology, i.e. what you are looking for and whether *Arabidopsis* is a good model of that investigation

Note, Springer Nature refused a link image request for the container overview, but the image I wanted to share is here.
> 2 votes

---
Tags: genome, software-recommendation
---


title: Check mutation of a gene

text:
# Question
Title: Check mutation of a gene
Hi all, I am checking if a gene which around 50k base pairs has mutation. So is that a variant from A to G? Would you please tell me how to do so? Thank you so much!

The reference sequence from GRCh38, sequenced using Illumina technology.

# Answer
Assuming this is a diploid genome, this looks like it could be a heterozygous variant (i.e. A is present on one chromosome, and G is present on the other chromosome).
> 1 votes

---
Tags: sam, vcf
---


title: Human genome similarities

text:
# Question
Title: Human genome similarities
As far as I know, we as humans share 99.9 similarity in our genomes. Further, we inherit 50% of the genome from our parents. So, it is correct to say that parent children have a 99.95% similarity in their genome?

# Answer
The answer is no if the question is read at face value,

1. I don't know the metric used (is this just SNPs?), nor which part of the genome is included in the calculation, just the exome?
2. Moreover, if the genetic identity between parents is unknown, then the similarity of their offspring to either parent cannot be calculated.
3. Your calculation excludes *de novo* mutations, particular at neutral sites and copy number variation, e.g. microsatellites, however its not clear whether the data is simply the exome.
4. Homologous recombination is excluded, i.e. particularly occurring in the parental germ-line, there must be stochastic effects even in the first generation.
5. Finally, mtDNA is maternally inherited.
> 2 votes

# Answer
The common wording is that the difference between two human *haplotypes* is 0.1%. This number only counts base substitutions, not insertions/deletions (indels) or structural variations (SVs). SVs together affects more bases than base substitutions. Also, the number is smaller is non-Africans (~0.07%) and higher in Africans (~0.12%). Anyway, let's denote this number by $\theta$, also known as the scaled mutation rate.

Because human is diploid, you need to count genotype differences when you compare two human individuals. Under the Wright-Fisher model, the rate of genotype difference on autosomes is

$$\theta\cdot\left(1+\frac{1}{2}\cdot\frac{1}{2}+\frac{1}{3}\right)=1.58\theta$$

The rate between a child and his/her parent is, I think:

$$\theta\cdot\left\[1\cdot\frac{1}{2}+\frac{1}{2}\cdot\left(\frac{1}{3}+\frac{2}{3}\cdot\frac{1}{2}\right)+\frac{1}{3}\cdot\frac{1}{2}\right\]=\theta$$

Note that I am sure the rate can be theoretically calculated but I am not confident in my calculation. Also, the rate of genotype difference on sex chromosomes or mitochondrial DNA is different.
> 2 votes

---
Tags: genome
---


# Question
Title: How do I calculate inbreeding on a single VCF file?
I would like to know the level of inbreeding for a single sample VCF file. I'm agnostic to the specific method of measuring inbreeding; any should work. However, I'd like to avoid having to specify a "background" or "reference" population. Ideally, the only input should be the single VCF file containing one human genome sample (WGS, hg38). (I will run this multiple times on different samples, but I'd like to run them individually.)

The existing tools seem to all require multiple VCFs.

Relatedly, is runs of homozygosity an acceptable metric for evaluating inbreeding in this context? And are there any special considerations to consider for various ethnic groups, like Ashkenazi Jews, that might influence the interpretation?

# Answer
**Possible solution**:

```
for record in vcf_reader:
    if not record.is_snp or len(record.ALT) != 1:
        continue  # skip if not a biallelic SNP

    gt = record.genotype('Sample')['GT']

    if gt == '0/1' or gt == '1/0':
        obs_hets += 1
    elif gt == '0/0' or gt == '1/1':
        obs_hom += 1

num_variants = obs_hom + obs_hets
exp_hets = num_variants * 0.5  # under the assumption of Hardy-Weinberg Equilibrium

f = 1 - (obs_hets / exp_hets)  # calculate inbreeding coefficient

```

Example output: "-0.25"

Interpretation: We observed more heterozygosity than we expected.

I don't know if this method is valid, however.

**A different way to calculate it**:

```
for record in vcf_reader:
    if not record.is_snp or len(record.ALT) != 1:
        continue  # skip if not a biallelic SNP

    gt = record.genotype('Sample')['GT']

    if gt == '0/1' or gt == '1/0':
        obs_hets += 1
        allele_counts['A'] += 1
        allele_counts['a'] += 1
    elif gt == '0/0':
        allele_counts['A'] += 2
    elif gt == '1/1':
        allele_counts['a'] += 2

total_alleles = sum(allele_counts.values())

# calculate allele frequencies
p = allele_counts['A'] / total_alleles
q = allele_counts['a'] / total_alleles

# calculate expected heterozygosity
exp_hets = 2 * p * q * total_alleles / 2  # the total number of variants

f = 1 - (obs_hets / exp_hets)  # calculate inbreeding coefficient

```

Example output: "-0.45"

Interpretation: We observed more heterozygosity than we expected. I don't have a good idea of the meaning of this number.
> 0 votes

---
Tags: vcf, genome
---


# Question
Title: How to determine the phylogeny of a specific set of genes
My goal is to understand the phylogeny of genes involved in bacterial restriction-modification systems, in psychrophilic bacteria.

To do so here is what I have done.

1. Downloaded the protein sequences of all RM genes in psychrophilic bacteria from REBASE. I also downloaded a set of genomes I'd like understand the phylogeny of from NCBI.
2. Clustered these genes using usearch's `uclust` method.
3. Aligned the resulting centroids using MUSCLE5.
4. Built an HMM from the aligned centroids using HMMER3.
5. Ran `hmmsearch` to get hits from my HMM on each genome (AA).

I am unsure how to proceed now. I now a phylogenetic tree requires an alignment, but I'm not sure where my HMM hits fit-in. Am I meant to use `hmmalign` on the output of `hmmsearch` and then feed the resulting alignment to IQ-TREE or similar?

Thank you for your help as I learn bioinformatics.

PS: Reading through the forum it seems like MOLE-BLAST does exactly this, however I have too many genes and their length is too long to use the MOLE-BLAST as hosted by NCBI at least.

I should mention I set the identity threshold of search to 50% for proteins. Also, I got this method from this paper which achieved something similar (https://pubmed.ncbi.nlm.nih.gov/30893937/).

# Answer
> 0 votes
Your strategy is fine but the order is wrong.

What I would do is:

* 2 -\> checked by 4 and 6

Run *usearch* and check the output of that in HMMER3. You can then use this to kick out stuff from usearch, that are not real homologues. Thats what HMMER3, hmmeralign (I dunno hmmersearch), are doing (below). It's also possible the HMMER approaches could replace usearch.

Once you have identified your homologues with confidence you can then align each separate group via *muscle5*

**Your way**

What you are doing is using HMMER3 to assess the quality of the alignment. What could be a problem is feeding *usearch* straight into muscle5 ... a non-homologous protein will get you a bad alignment. The bad alignment will not give you a good tree. *Both aligners and trees need nice homologues, otherwise they don't work.* The errors between aligner (muscle5) and tree are also synergistic. Thus if that bad sequences isn't removed by hmmer stuff then alignment errors compound onto sequence errors in the bad sequence giving a really bad tree.

You could of course align, HMMER3, filter stuff then align again (thats much better a bit complicated).

I don't know how good usearch is - but you are relying on it alot, supposing it gives you orthologues at \<25% identity that could be any protein. Basically you are feeding that stuff straight into muscle5 and if HMMER3 doesn't trap it that will go straight onto a tree. That could be a bad alignment and that would be a bad tree.

**Summary**

1. What you have to understand is phylogenetic trees are very sensitive to bad sequences, i.e. stuff that is not true homologues. The aligner will also be sensitive. There will also be bad synergistic relationships between the two.
2. HMMER3 and is associated methods (hmmeralign) are trying to clean stuff up before that stage, rather than simply relying on *usearch*. Prof Eddie (HMMER and pfam\[?\] author) has an excellent reputation in this area - not in tree building but identifying protein families and thats what you are doing here.
3. hmmeralign doesn't replace muscle5, but it will provide a very good assessment of the protein families around a given seed. It's going to give you the homologues. That output would then go to muscle5 to get homologous sites.
4. The tree needs *both* homologues and homologous sites.

The whole hmmer project has a really good reputation in the area of protein family identification and thats why its been recommended.

**Probable simple answer** I think the answer to the question is you're are confusing hmmeralign with muscle5. They are complimentary but distinct. Again hmmeralign identifies homologues and the muscle5 designates homologous sites.

---

Answering the comments

> 50% cut-off

It looks like *usearch* clusters according to a %age identity threshold.

There are three weaknesses:

1. Gene duplication can take place about 50% ... it depends whether thats cool, usually it isn't, but it's your call.
2. There can be complete change of function \>50% identity and this is about the definition of a homologue. If you are using a functional definition its not a homologue
3. Homologues can have \<\<50% identity and thats very common certainly in viruses I don't know about bacteria.

hmmersearch is going to rank your data according to the hmmer algorithms. It's going to be like blast,but all hmmer algorithms are focused around defining protein families (it's a better approach than blast *for this* application). So it provides another option to assess the output of *usearch*. Thus if there a gene duplication for example you should be able to spot it (it will be a sharp break in the continuous output). This will help you assess you 50% cut-off.

hmmeralign would be quite cool, because you can take a representation of *usearch* and then search that via hmmeralign and see if you get a similar output. It is extremely important for assessing 50% cut-off, does the homology drop below 50%? hmmeralign will tell you, but it must see the complete original database (not the stuff extracted from *usearch*. You can cross check that with hmmersearch to assess where hmmeralign is making the cut-off. So hmmeralign will independently (or *usearch*) will find you homologues and hmmersearch will give you an insight about how it's done it.

---
Tags: phylogenetics, phylogeny, hmmer
---


# Question
Title: Pal2nal translation of large multi-fasta files produces a codon translated file where some sequences half length of the average
I did sequence alignment of a large peptide multi-fasta (n= 4991 sequences). The peptide alignment has sequences with the same length and pal2nal went through just fine... except some of the codon sequences are at half length. If average is X then some sequences are X/2. This is choking IQ-Tree.

I have tried both MUSCLE super5 and MAFFT. The error remain the same (i.e. MUSCLE or MAFFT both lead to some sequences having half of average length) except for different average lengths and average half length in MUSCLE and MAFFT codon sequences. I have pulled out and played with the sequences causing the issue and they seem to be in frame.

Example of peptide sequence not causing an issue: `RKVEAFLLFKEMGERGCQPNVHTYTVLIDSFCKERNLDDARKLFDDMFKKGLVPSVVTYNALIDGYCKEGMTEAALEILGMMESKKCNPNARTYNELICGFCKAK`

corresponding cds `AGGAAAGTGGAAGCTTTTCTACTTTTTAAAGAAATGGGTGAAAGAGGTTGTCAGCCTAATGTTCATACATACACTGTGCTTATTGATTCCTTCTGTAAGGAAAGGAATCTTGATGATGCCAGGAAATTGTTTGATGACATGTTTAAGAAAGGTTTGGTTCCCAGTGTGGTCACTTATAATGCTTTAATTGATGGGTATTGTAAAGAGGGAATGACTGAAGCTGCATTAGAAATTTTAGGTATGATGGAATCAAAGAAATGCAACCCTAATGCTCGGACCTACAATGAATTGATCTGTGGATTTTGTAAAGCTAAA`'

Example of peptide causing issue:

`GLCKGGRLNDAWEIFQYLLAKGYQLNVHTYNAMVHGFCKEGLLDEAISLLYKMEENGCVPNSVTFNVVL`

corresponding cds `GGTTTGTGCAAAGGTGGTAGATTAAATGATGCGTGGGAGATTTTTCAGTATCTTTTAGCGAAAGGTTATCAACTAAATGTCCATACATATAATGCGATGGTTCATGGTTTTTGCAAAGAAGGTTTGCTTGATGAAGCAATCTCCCTGCTTTATAAAATGGAAGAGAATGGTTGTGTCCCTAATTCTGTAACTTTTAATGTAGTCCTT`

Any idea what might be going on?

# Answer
> 1 votes
I had a quick blast of the example sequences and this is the *Acacia* so ancient vascular plant. It forms a locus within the protein in the first part of the question. This protein is comprises a large repeat amino acid motif (check NCBI), the second "issue sequence" described in the question was close to the size of the repeat locus, but start/stop of the locus is not in accordance with the repeat structure. Thus I don't think the repeat locus is part of the issue.

**Genome assemblies** I noted the examples you gave were all draft genomes and this suggest the problem was with the original genome assemblies. If this was true the way to overcome this would be to reassemble the sequence using a reference genome, *bwa-mem2* type thing. Thus, here the error is occurring way before the sequences is even downloaded. Just to note, repeat structures can often be collapsed on a draft genome assembly, but I don't think thats happened here.

**Reading frame** what I think you are doing is performing an alignment in the proteins and then using that alignment as the template which is then transferred to the sequence data. So now the nucleotide data is nicely aligned including all indels. That appears to be the goal of *pal2nal*. A common error here is reading frame error, particularly at the beginning of the sequence. I use *transalign* regularly for this purpose and do encounter that error. I don't know how *pal2nal* deals with reading frame errors so I can't comment, but thats why I wanted to see the nucleotide input to *pal2nal* to check this.

**Trees** the reason this is taking place is because *iqtree* doesn't permit a tree. That I can help with, I know about tree building. If you post your *iqtree* command as a separate question , together with the *iqtree* error, I can probably answer that. What I'll identify is how to build a tree, this likely to exclude sequence, or could be addressed as two trees. Dunno, but if the bottleneck is can't build a tree - thats easily fixable. In this case, is the sequence issue an assembly error, reading frame error ... something else. There's not enough data that I could point to as the underlying source.

---
Tags: phylogenetics, sequence-alignment, phylogeny, proteins, codon
---


# Question
Title: Keeping DNA sequence after changing FASTA header on command line
I have a FASTA header that looks like this:

```
>7c8250ef-c89f-4d42-9d48-12c8fe245fb2 runid=606f271fc97598006ba5a922136a2c304cef75a5 sampleid=Pool12-1 read=19008 ch=301 start_time=2021-07-03T08:48:18Z barcode=barcode01

```

And I am able to change it to the desired output here:

```
>7c8250ef-c89f-4d42-9d48-12c8fe245fb2_001

```

Using this command:

```
grep '^>' 001_old.fasta | cut -d ' ' -f 1,8 | sed '/^>/s/$/_001/' > 001_new.fasta

```

However it completely deletes the DNA sequence for that read and only outputs the names of the headers to the new fasta file:

```
>7c8250ef-c89f-4d42-9d48-12c8fe245fb2_001
>815440d5-e033-446c-b7dd-5fad4105184a_001
>d36a0fea-101d-4bd9-bee4-adec7fe4d43c_001
>38ca9eb8-bbd2-4d0a-ab06-79d134ad46f1_001
...

```

I think it has something to do with the `cut` command I'm using. I tried just doing 1 but that gave the same output so I figured if I included the 8th position of text after a space it would work but again no dice. If someone could help me to include the DNA sequence under the FASTA header in this new file that would be great.

# Answer
> 5 votes
The `grep`, `cut`, `sed` pipeline could also be replaced with the following command:

```
awk '/^>/ { print $1 "_001"; next }1' file.fa

```

The general syntax of an AWK command is `condition { action }` where `condition` specifies a condition or regular expression, and `action` specifies a block of code or set of commands to be run when the condition is met. Both `condition` and `action` are *optional*. Here, we use a regular expression to match lines beginning with a "\>" symbol, and print the first whitespace-separated field followed by "\_001". We then use `next` to skip the remaining actions and move to the next line, which avoids any further processing of the current line. The `1` on the end is actually another condition that always evaluates to true. Since there's no action associated with this condition, AWK performs the default action which is to print the entire line. This effectively prints all lines that do not begin with a "\>" symbol.

# Answer
> 3 votes
The problem is in the `grep` command you're using, which only sends the fasta file headers to the output. Here's my guess at what the commands in your pipeline do:

1. `grep '^>' 001_old.fasta` \- extract lines starting with '\>' from `001_old.fasta`
2. `cut -d ' ' -f 1,8` \- assume file contains space-delimited records; extract only the first and eighth field from each line
3. `sed '/^>/s/$/_001/'` \- in lines that begin with a '\>', replace the end of the line with '\_001'
4. `> 001_new.fasta` \- write result out to `001_new.fasta`

Of all of these, I think only the `sed` and output parts are important, but it needs to be slightly modified to work as intended:

```
sed '/^>/s/ .*$/_001/' 001_old.fasta > 001_new.fasta

```

Assuming I've got the syntax right, this should search for lines beginning with '\>', and replace everything from the first space up to the end of the line with '\_001'.

I'm more familiar with perl, so here's a perl one-liner that should do the same thing:

```
perl -pe '/^>/ && s/ .*$/_001/' 001_old.fasta > 001_new.fasta

```

---
Tags: fasta, grep, multi-fasta, command-line
---


# Question
Title: Correlation heatmap of RNA-seq clusters all samples together leading to very low no. of DEGs
I am writing to you to take an input or may be you can provide a different perspective. I am at wits end :(

So I have nearly 200 samples. These are separated into two groups (group I treated with drug1 & group II treated with drug II) of course these two groups have sub groups of responders to drug and non responders.

I have used DESeq2 to analyze this data. The first thing that I notice about this data is that the in the correlation plot the samples clustering across the groups. In past I have seen some samples clustering across the groups but this a lot. There is no cluster specific to one group.

Secondly, I get like 5 or 6 DEG that are statistically significant and have fold change of 2.

I have been part of similar studies with same tissue and comparable drugs but have never noticed something like this.

I am not sure what is happening here and why it is happening.

Any help or insights are sincerely appreciated.

# Answer
> 4 votes
This answer has two stages,

1. Diagnostics, what analytically has gone wrong?
2. Where to go from here.

I am pretty confident about point 1. I've been in comparable situations and I know trees. Stage 2, is more open to opinion, however there's no point criticising something (point 1) without pointing out how it can be corrected (point 2).

If your question is within group heterogeneity, excepting for the large group of non-responders associated with D1, thats a different answer.

Low number of DEGs is about your control group? Difficult to comment because the DEG calculation is not declared, nor volcano plots, cut-offs, FDR etc...

---

You have to keep in mind this data mining and the there are some interesting groups emerging. For first stage analytics looks reasonable (low DEGs dunno).

**Result** You have 9-10 clusters and those are informative. You can see this because they are forming a checker-board across the heat-map.

**Issue** The issue is that the deeper branches below the 9-10 clusters are not informative. Do see? The deep branches are likely to be repositioned in different orders in different analyses. Those alternative order might make better biological sense. On the heat-map shown there's a risk of attaching two much important to deep branching orders for a UPMGA tree for a correlation dendrogram.

* Heterogeneity across the matrix is not permitted in UPGMA and causes artefacts - because you are trying to say there's a uniform rate when in most biology that assumption doesn't work. This will usually result in incorrect branching orders.
* You are forcing the data into a bifurcating tree: it might not be at all, could be a large polytomy (star tree).
* Even data sets designed for trees (correlation matrixes are not good data sets for trees) ... branching orders go wrong for loads of reasons.

The second stage analysis would be to take the correlation matrix and apply the neighbour-joining (nj) algorithm to it and see whether the deeper branch align differently. There are also *much* better tree methods, but that starts getting complicated fast, because there are alternative approaches, such as co-variance. Advanced tree methods (Bayes and ML) will model shifts heterogeneity - however it needs some expertise to get it working.

I'd start assessing the assumptions UPGMA vs neighbour-joining does that take me closer to a biological explanation, then a 9x9 checkboard? Thus would the number of squares in a checkerboard by greatly reduced via nj?

**Machine learning** Supervised learning (ML), is the trendy way into the data set and in this day an age thats probably the way I'd go. What I'd do is use the 9-10 clusters here as an initial assessment of the training targets. I'd then perform whats called a "feature selection analysis" and assess whether that output makes biological sense. The feature selection should start identifying what parts of your data are causing the shifts between those 9-10 clusters.

---
Tags: rna-seq, phylogenetics, phylogeny, deseq2
---


# Question
Title: Seurat: Error in FetchData.Seurat(object = object, vars = unique(x = expr.char[vars.use]), : None of the requested variables were found:
I was trying to do subsetting via:

```
pbmc = subset(pbmc, subset = nFeature > 1000 & nCount > 4000)

```

R returns:

```
Error in FetchData.Seurat(object = object, vars = unique(x = expr.char[vars.use]), :
None of the requested variables were found:

```

and then nothing, so it says cannot find requested variables and doesn't tell me what's missing...

I'm very confused. Thank you so much for your help!

# Answer
I've come across the same error which seems to be a common issue when using the `subset()` function with Seurat objects. Sometimes, I encountered the same error and couldn't figure out the reason behind it, checking the variable's names was not a remedy in my case as I specified the exact same name for subsetting as to what I could see in my object metadata. However, I found a workaround to subset cells in a Seurat object using a different approach.

Instead of using `subset()`, you can try the following:

```
seurat_obj_subset <- seurat_obj[, <condition to be met>]

```

For example, if you want to subset a Seurat object called 'pbmc' based on conditions like having more than 1000 features and more than 4000 counts, you can use the following code:

```
pbmc = pbmc[, pbmc$nFeature > 1000 & pbmc$nCount > 4000]

```
> 1 votes

# Answer
From the steps in the PBMC\_3k tutorial, it looks like you need to suffix the name with the assay as well, separated by an underscore (in this case, `RNA`):

```
pbmc <- subset(pbmc, subset = (nFeature_RNA > 1000) & (nCount_RNA > 4000))

```
> 0 votes

# Answer
> "Error in FetchData.Seurat(object = object, vars = unique(x = expr.char\[vars.use\]), : None of the requested variables were found:"

As the error is pointing out here, you should check your variable names that you have used here. It must match the variable names in the pbmc object that you want to subset.
> 0 votes

---
Tags: r, seurat, subset
---


# Question
Title: How can I give different names to files in a directory with a for loop in a bash script?
I'm expecting to get 17 different paired-end fastq files (34 in total), so I want to make a bash script to just run my code through all the fastq files in a directory at once. How can I change the name of the input and output files each time the script runs through each file? So when it moves to the file\_002, all names have file\_002 at the beginning instead of file\_001, and so on. And also, when merging the R1 and R2 reads how can I make that it only merges the correspondant files with a loop? for examples merging only file\_001\_R1 with file\_001\_R2, file\_002\_R1 with file\_002\_R2, file\_003\_R1 with file\_003\_R2, and so on.

This is my code:

```
for file in directory_name
do
pear -f file_001_R1.fastq.gz -r file_001_R2.fastq.gz -o file_001.fastq
cutadapt -g TGATAACAATTGGAGCAGCCTC...GGATCGACCAAGAACCAGCA -o file_001_barcode.fastq file_001.fastq
cutadapt -g GTGTACAAATAATTGTCAAC...CTGTCTCTTATACACATCTC -o file_001_UMI.fastq file_001.fastq
seqkit concat file_001_barcode.fastq file_001_UMI.fastq > file_001_concatenation.fastq
seqkit rmdup -s file_001_concatenation.fastq -o file_001_unique_pairs.fastq
seqkit subseq -r file_001_unique_pairs.fastq > file_001_unique_barcodes.fasta
bowtie -q --suppress 1,2,4,6,7,8 -x ref_index file_001_unique_barcodes.fasta > file_001_barcodes_allignment.bowtie
sort file_001_barcodes_allignment.bowtie | uniq -c > file_001_barcode_counts.txt
awk 'BEGIN{print "Barcode,TF_variant,Code"}{print $3","$2","$1}' file_001_barcode_counts.txt > file_001_barcode_counts.csv
done

```

Thank you.

# Answer
Here's what I use - along with an example you can modify as needed:

```
FASTQ_DIR="/directory/containing/FASTQs/"

for R1_FASTQ in "$FASTQ_DIR"*R1*.fastq.gz
do
                         BASE=$(basename "$R1_FASTQ" _R1.fastq.gz)
                         R2_FASTQ="${R1_FASTQ/R1/R2}"
                         R3_FASTQ="${R1_FASTQ/R1/R3}"
                         echo "BASE: $BASE"
                         echo "R1: $R1_FASTQ"
                         echo "R2: $R2_FASTQ"
                         echo "R3: $R3_FASTQ"
                         # process R1/R2 here using "$R1_FASTQ" "$R2_FASTQ" for input, and "${BASE}_XXXXXXX.fastq" for output
done

```

Basically, the strategy is to iterate over all the R1 FASTQs in a directory, calculate a `BASE` name that can be used to rename output files with arbitrary extensions (like `.sam`, `.bam`, etc...), and use BASH shell parameter expansion to calculate and set variables to the R2 (and R3 if necessary) reads present in an Illumina run that uses standard Illumina naming conventions.

This may need modification depending on the naming conventions you are using.
> 3 votes

# Answer
Something like this? This uses `printf` to create a zero-padded 3-digit number corresponding to 17 numbers from 1 to 17:

```
for fileBase in $(printf 'file_%03d\n' $(seq 1 17));
   do echo ${fileBase}_R1.fq.gz ${fileBase}_R2.fq.gz;
   pear -f ${fileBase}_R1.fastq.gz -r ${fileBase}_R2.fastq.gz -o ${fileBase}.fastq
   ...
done

```
> 2 votes

---
Tags: bash
---


title: Loop Pymol command get_area

text:
# Question
Title: Loop Pymol command get_area
I want to select multiple residues from different PDB files loaded on the same pymol session.

```
select resi LIG
get_area sele

```

How do I loop the above command to enable the command below:

```
select ligands,resn LIG

```

# Answer
> 3 votes
from https://www.blopig.com/blog/2019/08/how-to-iterate-in-pymol/

Example 1: Print the names of all the selected objects

You’ve selected several objects in PyMOL, and you want to see their names. By default, a new selection object called “(sele)” is created. We can create a Python variable called seleobjs, say, and then loop over them:

```
seleobjs = cmd.get_object_list('(sele)')
    for obj in seleobjs: print obj

```

Example 2: Disable all the selected objects

If you still have your selection, “(sele)”, you can disable them (effectively turn them off on the right-hand sidebar):

```
seleobjs = cmd.get_object_list('(sele)')
    for obj in seleobjs: cmd.disable(obj)

```

Example 3: Delete all the selected objects

You’re done with the objects in your selection, “(sele)”: you might be able to delete them using “delete ”, but if you can’t write a wildcard that covers all the names of the objects, you can iterate over a selection and delete them:

```
seleobjs = cmd.get_object_list('(sele)')
    for obj in seleobjs: cmd.delete(obj)

```

from https://pymolwiki.org/index.php/Get\_Area

get\_area calculates the surface area in square Angstroms of the selection given. Note that the accessibility is assessed in the context of the object(s) that the selection is part of. So, to get the surface areas of e.g. a component of a complex, you should make a new object containing a copy of just that component and calculate its area.

please let us know if it worked

# Answer
> 1 votes
Ok went back to this question having `pymol` within a hand's reach this time. Question is not really clear to me so I am posting a reasonable (to me) version of the question into a new code snippet;

I am using two PDBs containing `ATP` as ligand `'3dgl.pdb', '2p09.pdb'`

code:

```
import pymol

list_pdbs = ['3dgl.pdb', '2p09.pdb']

pymol.finish_launching()

for _i in list_pdbs :

    pymol.cmd.load( _i, _i.rstrip('.pdb'))

objs = pymol.cmd.get_object_list(selection='(all)')

print('structure loaded : ' , objs)

pymol.cmd.remove('resn HOH')

for _i in objs:

    pymol.cmd.create(_i+'_ATP_copy' ,  _i + ' and resn ATP')

    # print('after create ')

    area = pymol.cmd.get_area(_i+'_ATP_copy', quiet = 0)

    pymol.cmd.show('dots' , _i+'_ATP_copy')

    print(_i+'ATP_copy' + ' area --->  ' + str(area))

```

output:

```
structure loaded :  ['3dgl', '2p09']
 cmd.get_area: 1013.704 Angstroms^2.
3dglATP_copy area --->  1013.7040405273438
 cmd.get_area: 1013.704 Angstroms^2.
2p09ATP_copy area --->  1013.7040405273438

```

picture of the `ATP` surfaces in Pymol itself,

you can play with `cmd.get_area(...)` parameters see https://pymolwiki.org/index.php/Get\_area

---
Tags: python, pymol
---


# Question
Title: How do I add reference/alternate alleles to my vcf file?
Here are the first few columns of my VCF v4.2 file. As you can see, the alt allele is missing. The VCF's refs and alts don't correspond to the reference genome.

| #CHROM | POS | ID | REF | ALT |
| --- | --- | --- | --- | --- |
| 9 | 77051 | rs9408670 | T | C |
| 9 | 114350 | rs12685480 | A | . |
| 9 | 116581 | rs558835286 | G | . |
| 9 | 117550 | rs267602077 | A | . |
| 9 | 117720 | rs10959293 | T | . |
| 9 | 126903 | rs62533602 | G | A |
| 9 | 133828 | rs78804900 | C | T |
| 9 | 154736 | rs200683662 | T | . |
| 9 | 154933 | rs17795008 | T | G |
| 9 | 174425 | rs2998944 | T | C |
| 9 | 174500 | rs142581442 | G | . |
| 9 | 175050 | rs62532020 | G | A |
| 9 | 185632 | rs2992854 | C | A |
| 9 | 188894 | rs3951559 | G | A |
| 9 | 191843 | rs3008128 | G | A |
| 9 | 202061 | rs117675283 | T | . |
| 9 | 204201 | rs10964134 | G | A |

The missing alt allele happens when the genotype is 0/0. The VCF file isn't saying anything incorrect. However, there is a tool I'm trying to use that filters out "0/."genotypes.

A better way of formatting them (which I'm trying to achieve) would have the ref/alt alleles correspond to the reference genome, and just have a redundant letter in the case of a homozygous genotype.

I'd like my VCF to display the correct alt/ref for that position even if there are no samples with a genotype corresponding to the alt allele. How do I do this? I'm not sure what this process is called.

Ideal output:

| #CHROM | POS | ID | REF | ALT |
| --- | --- | --- | --- | --- |
| 9 | 77051 | rs9408670 | T | C |
| 9 | 114350 | rs12685480 | A | G |
| 9 | 116581 | rs558835286 | G | A |

# Answer
> 1 votes
Let me know if you see any potential issues with this solution. It just replaces the dots with random letters (easier/faster than aligning to reference genome).

```
def add_ref():
    # Directory path
    dir_path = '/path/'

    # Function to process genotypes in each row
    def process_genotypes(row):
        genotypes = np.array([x.split(":")[0] for x in row[9:]])
        if row[4] == "." and not any(val in ["1/1", "0/1", "1/0"] for val in genotypes): # to make sure I'm not introducing an allele that isn't present
            # also can do if row[3] in ["A", "T"]: then random.choice(["C", "G"])
            if row[3] == "A":
                row[4] = "C"
            elif row[3] == "T":
                row[4] = "G"
            elif row[3] == "C":
                row[4] = random.choice(["T"]) # I know
            elif row[3] == "G":
                row[4] = random.choice(["A"])
        return row

    # Get all .vcf.gz files in the directory
    files = glob.glob(dir_path + '*.vcf.gz')

    for file in files:
        print(f"Processing {file}...")
        # Name of the new file
        new_file = os.path.join(dir_path, 'new_' + os.path.basename(file))
        with gzip.open(file, 'rt') as in_f, gzip.open(new_file, 'wt') as out_f:
            for line in in_f:
                if line.startswith('#'):
                    out_f.write(line)
                else:
                    row = line.strip().split('\t')
                    new_row = process_genotypes(row)
                    out_f.write('\t'.join(new_row) + '\n')

```

---
Tags: vcf, genome
---


# Question
Title: How do you extract the full list of files in a GEO archive?
I am new to use of GEO and having trouble figuring out a straightforward way to grab the names of the files, (through programmable handles like GEOparse or otherwise) that is the entire and complete list of all files fastq and otherwise (NOT downloading the actual data, but the names and descriptions) within a specific GEO archive (GSEXXX). Albeit with the intention of categorizing the data into more understandable and friendly to browsing segments, the data presentation in a GEO archive seems very fragmented, making it hard to find what you are looking for. When you submit your data to GEO, they require you to have a spreadsheet, that is essentially the list of files you are submitting with their corresponding meta-data descriptions. Is it possible to get that exact submission spreadsheet associated with a data archive?

# Answer
> 1 votes
Unlikely that you will be able to get the original submission spreadsheet, but you can use Geofetch

---
Tags: chip-seq, geoquery, geo
---


# Question
Title: Retrieving NCBI Taxa IDs from refseq or GenBank assembly accession
I have about 10,000 genome files all named by either refseq or genbank accession number, do you know if it's possible to convert these numbers to the corresponding NCBI taxon ID or species?

for example:

```
GCA_000005845.2 to 79781

```

In the case of *E.coli*

Edit:

The file names look like this:

```
GCF_000337855.1_ASM33785v1_protein.faa.gz
GCF_000091125.1_ASM9112v1_protein.faa.gz
GCF_000184535.1_ASM18453v1_protein.faa.gz

```

My operating system is Ubuntu 16.04

# Answer
Turns out I'd already written some code that did it, human memory is a funny thing.

This uses biopython to split the field description to where the species is. May not work for all NCBI files, but seems to work on most.

```
import Bio
from Bio import SeqIO
from Bio import AlignIO

for record in SeqIO.parse (FILE, "fasta"):
    Speciesname = record.description.split('[', 1)[1].split(']', 1)[0]

```
> 1 votes

# Answer
You can use eutils to get this information as shown below:

```
esearch -db assembly -q 'GCF_000337855.1' | esummary | xtract -pattern DocumentSummary -element AssemblyAccession,Taxid
GCF_000337855.1 469382

```
> 2 votes

# Answer
As the question mentions python, here is a python example using biopython.

```
import json
from Bio import Entrez

def accession2taxid(acc: str, db="nucleotide") -> str:
    handle = Entrez.esearch(db=db, term=acc)
    record = Entrez.read(handle)
    gi = record["IdList"][0]
    handle = Entrez.esummary(db=db, id=gi, retmode="json")
    result = json.load(handle)["result"]
    taxid = result[gi]["taxid"]
    return str(taxid)

```

So if you want to use an assembly accession, as in your example

```
>>> accession2taxid("GCA_000005845.2", db="assembly")
'511145'

```

In your example, you said that `GCA_000005845.2` should return `79781`, but that Taxon ID is for "uncultured archaeon AOC672". The taxid in the example I just showed returns "Escherichia coli str. K-12 substr. MG1655", which is what you want.

Your question also mentions GenBank accessions, so here is another example

```
>>> accession2taxid("CU458896.1", db="nucleotide")
'561007'

```
> 1 votes

---
Tags: python, phylogenetics
---


# Question
Title: What would be the best method to obtain every prokaryotic psychrophile genome?
I am interested in running an analysis on the genome of every psychrophilic bacteria. Here defined as bacteria that have an optimal growth temperature at less than 15°C.

What would be the best way to go about collecting these?

I have yet to find a "growth temperature" filter on databases like NCBI. Though REBASE does include that (however obtaining the genomes from that site does not seem straightforward).

Thank you

# Answer
> 2 votes
I ended up copying the GenBank (primary) column from the REBASE search results into an excel sheet. From there, I created a file with one accession per line called `Psych_Accessions.txt`. I then ran the following command using NCBI's command line tools.

```
epost -db nuccore -input Psych_Accessions.txt -format acc | efetch -format fasta > psych_genomes.fasta

```

# Answer
> 2 votes
The database to use is BactoTraits. This will help either in the entirety or at least provide a solid starting point. Its here

BactoTraits – A functional trait database to evaluate how natural and man-induced changes influence the assembly of bacterial communities by Aurélie Cébron (2021) Ecological Indicators 130

> We have created a user-friendly database of bacterial functional traits, thanks to the properties of **19,455 bacterial strains**. This database has been called BactoTraits. For example, oxygen preference, size and shape of bacteria, motility, optimum and range of pH and **temperature**, genome GC percent and trophic type are among the 19 traits included in BactoTraits. Based on the best-informed strains in the database, we identified five functional groups (i.e. mesophiles, competitors, colonizers, stress-sensitives and stress-tolerants) exhibiting a wide strain taxonomic diversity but with quite similar trait profile combinations.

Thus it's nearly 20000 species so your batteria should be there.

You can download the exact database here. This is Univ de Lorraine's site, France. Its 9.4MB and called

> BactoTraits\_databaseV2\_Jun2022.csv

The supplementary information in the paper gives all the details for interpreting the database columns.

Given your skills in development you should have no trouble converting the **csv** format into a functional database and interrogating/querying it. My understanding of the paper is that the authors are microbiologists rather than coders, which in context is perfect.

---

Linking NCBI

BactoTraits contains two strain designations

So ATCC links to NCBI links to the SRA and genome.

DSM are not necessarily available on NCBI, but ATCC will there. NCBI do have lots of databases for download, a database comprising ATCC to NCBI genome ... dunno, but *Entrez* will be able to link the ATCC (American Type Culture Collection)

---
Tags: metagenome, wgs, bacteria
---


title: Plotting infercnv results

text:
# Question
Title: Plotting infercnv results
I'm working with matched single cell data, where we have treated and untreated samples for the same patient. I ran CNV analysis using the `infercnv` package.

I've followed the tutorial:

```
# data matrix
counts_matrix <- scData@assays$RNA@counts
meta = data.frame(labels = Idents(scData), row.names = names(Idents(scData)))
unique(meta$labels) # check the cell labels

```

```
# Create the infercnv object
infercnv_obj <- CreateInfercnvObject(raw_counts_matrix = counts_matrix,
                                     annotations_file = meta,
                                     delim = "\t",
                                     gene_order_file = "hg38_gencode_v27.txt",
                                     ref_group_names = NULL)

```

```
# perform infercnv operations to reveal cnv signal
infercnv_obj = infercnv::run(infercnv_obj,
                             cutoff=1,  # use 1 for smart-seq, 0.1 for 10x-genomics
                             out_dir="output_dir",  # dir is auto-created for storing outputs
                             cluster_by_groups=T,   # cluster
                             denoise=T,
                             HMM=T
                             )

```

Now, I'm not sure how to plot the results. Any help would be appreciated.

When running the below recommendation I get this:

```
INFO [2023-05-19 12:49:22] ::plot_cnv:Start
INFO [2023-05-19 12:49:22] ::plot_cnv:Current data dimensions (r,c)=34108,17624 Total=82463956 Min=0 Max=55099.
INFO [2023-05-19 12:49:22] ::plot_cnv:Depending on the size of the matrix this may take a moment.
<sparse>[ <logic> ]: .M.sub.i.logical() maybe inefficient
Warning: sparse->dense coercion: allocating vector of size 4.5 GiBINFO [2023-05-19 12:50:12] plot_cnv(): auto thresholding at: (-2.725632 , 3.000000)
INFO [2023-05-19 12:50:35] plot_cnv_observation:Start
INFO [2023-05-19 12:50:35] Observation data size: Cells= 17624 Genes= 34108
INFO [2023-05-19 12:50:35] clustering observations via method: ward.D
INFO [2023-05-19 12:50:35] Number of cells in group(1) is 331
INFO [2023-05-19 12:50:36] group size being clustered:  331,34108
Error in parallelDist(data_to_cluster, threads = infercnv.env$GLOBAL_NUM_THREADS) :
  x must be a matrix or a list of matrices.

```

# Answer
Personally I would not use `intercnv_obj` throughout, because is there was a bug in line 2, i.e. using `=` instead of `<-` then the created `intercnv_obj` might be trying to be passed into `plot_cnv`.

Could you try

```
data(infercnv_obj)

plot_cnv(
  infercnv_obj,
  out_dir = ".",
  title = "inferCNV",
  obs_title = "Observations (Cells)",
  ref_title = "References (Cells)",
  cluster_by_groups = TRUE,
  cluster_references = TRUE,
  plot_chr_scale = FALSE,
  chr_lengths = NULL,
  k_obs_groups = 1,
  contig_cex = 1,
  x.center = mean(infercnv_obj@expr.data),
  x.range = "auto",
  hclust_method = "ward.D",
  custom_color_pal = NULL,
  color_safe_pal = FALSE,
  output_filename = "infercnv",
  output_format = "png",
  png_res = 300,
  dynamic_resize = 0,
  ref_contig = NULL,
  write_expr_matrix = FALSE,
  write_phylo = FALSE,
  useRaster = TRUE
)

```

---

If the above doesn't work ...

What I think is happening is the input data is the problem. This precise error is already known and was reported on Broad's GitHub here:

https://github.com/broadinstitute/infercnv/issues/415

The responder (Broad Institute) believed it was the version problem, i.e. the version was outdated. However a separate respondee believed it was the input data format and suggested shifting away from a sparse matrix. That would make sense in context to the error being thrown.
> 0 votes

# Answer
If you haven't already, I recommend checking out the infercnv YouTube tutorial (https://github.com/broadinstitute/inferCNV/wiki/Running-InferCNV) for detailed instructions on running InferCNV. Additionally, you can explore Uphyloplot2 (https://github.com/harbourlab/uphyloplot2) for building a tree from the output generated by InferCNV.

I noticed that you didn't use the untreated sample as the normal reference for your CNV analysis. I'm curious to understand your reasoning behind this. In my case, I also have an untreated sample that isn't diploid, along with two treated samples. Initially, I planned to use the untreated sample as the "normal" reference. However, considering that InferCNV assumes a normal diploid reference, I'm now uncertain if using the untreated control would be appropriate.

Please let me know your thoughts, any help would be much appreciated!
> 0 votes

---
Tags: r, single-cell, cnv
---


# Question
Title: Statistical methods suitable for DE analysis of non coding RNAs
I am currently working on DE analysis of coding as well as non-coding features of *A. thaliana* using the `edgeR` package.

Is the negative binomial method that is normally used for the DE analysis of genes suitable for non-coding features such as miRNA, lncRNA, snRNA, snoRNA, pseudogenes, transposons, uORFs and 5' and 3' UTRs? Or are there special statistical methods that are more suitable for modelling each of these, like I have read that the zero-inflated model used by `lncDIFF` can better model lncRNA counts? I see that my miRNA counts matrix is abundant in zero counts and very low counts of genes. Should I use lncDIFF for that?

Also, among the various DE test methods in edgeR such as `exactTest()`, `glmFit()` and `glmLRT()`, `glmQLFit()` and `glmQLFTest()`, which is/are the one(s) most suitable for running DE tests on such non-coding features?

# Answer
> 1 votes
There is no question that **lncDIFF** is the approach to use, especially if the alternative is performing **edgeR** directly, and thats demonstrable (below). I don't think anyone would question the legitimacy of your approach. If you look at Fig 1 for the paper describing the test, here there is no question about it **edgeR** used directly does not perform for this task at all. The paper describing lncDIFF is here,

https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-019-5926-4

You present **lncDIFF** as an option within **edgeR**. This I don't know, but **lncDIFF** can be directly downloaded from GitHub, https://github.com/qianli10000/lncDIFF.

---

*The one thing you have to be careful about is you are using plants, whilst these tests are designed around humans, see below for details.*

---

How it's tested: GLM methods versus an Exact test, e.g. Fisher's. Both are all valid tests. However, GLM methods are preferred over an exact test, in my statistical experience.

First, you need to obtain a `fit` before you can perform a test, which ever GLM you choose.

```
x <- model.matrix(...)
y <- estimateDisp(...)
fit <- glmFit(y, x, robust=TRUE)
lrt <- glmLRT(fit, contrast=c(-1,1,0,0,0))

```

The basic issue is GLM likelihood ratio test versus quasi likelihood test, which the right approach? The clear answer is QLT because the authors of **lncDIFF** developed the test around QLT (see the paper). That shouldn't be general advice outside **lncDIFF**

The correct approach is to try both, if there major differences thats an issue: I would say in that scenario is consider QLT - but only for **lncDIFF**, as described. If there are no major differences you are good to go.

---
Tags: rna-seq, differential-expression, edger, non-coding
---


# Question
Title: Filtering criteria for non-coding features with very low counts
I am trying to do DE analysis of non-coding features of *A. thaliana*. I find in the miRNA and lncRNA counts file that they are abundant in zero counts, and most of the non-zero counts are very low. Now, for the DE analysis of coding features, an inbuilt filtering method to remove the extremely low counts is available in edgeR, called `filterByExpr()`. Is this function suitable for only those features with high counts like coding sequences, or is it suitable for non-coding features with low counts and low library sizes too? If not, then what kind of filtering criteria should I use for filtering out genes with insignificant counts? Or should I filter out only the ones with row sums 0?

# Answer
> 0 votes
Partially answered before. I am assuming this is a **lncDIFF** calculation. For the **lncDIFF** paper it is unlikely to work as the following figure displays,

https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-019-5926-4/figures/4

However, don't take my word for this, test it. You should test the option explicitly via `glmQLFTest` and `glmLRT`. If you perform either of these tests with or without the filter the better fit (one not rejected) is the correct answer. I strongly predict the filtered data will be rejected and unfiltered data (hopefully) will pass.

The question of what is better `glmQLFTest` versus `glmLRT` is discussed in the previous answer, here Statistical methods suitable for DE analysis of non coding RNAs

Just keep in mind this is plant data and these tests are developed for humans.

---
Tags: rna-seq, differential-expression, filtering, edger, non-coding
---


# Question
Title: is it possible to count cell in Violin plot in seurat?
I tried to quantify the number of cells of Fzd9+ Macrophage.

how can I know the cell count?

thank you

# Answer
> 3 votes
This can be a good starting point: `table(your_seurat_object@meta.data$the_column_that_has_the_labels_of_interest)`

Basically you will be counting the number of rows (each single cell is a row) in the `meta.data` slot of your Seurat object for each level of your column of interest (for example cell type annotations).

---
Tags: rna-seq, scrnaseq
---


# Question
Title: Issue with Differential Gene expression analysis with Deseq2
I am using the Deseq2 package to perform a DE analysis on multiple factors.

So, basically to perform the analysis I consider patients characterized by the Brugada Syndrome type 1 and other patients that instead are controls.

In order to have equal numbers of PZ and CTRs I included in my analysis all the PZ of type 1, that in this case are 8 samples, and I selected 8 random CTRs among the full dataset the CTRs that I have for this batch (in total are 40).

To select the random CTRs I used the sample\_n function:

```
random_CTR_samples_IGATECH<-sample_n(samples_IGATECH_ALL_CTR,8)

```

Anyway I created the design matrix that looks like this including also the Age variable as well as the Gender one since I want to test also the effect of those in my data, (NON\_INDUCIBLE)= CTR:

```
       Sample BrS_BASELINE_PATTERN Gender Age
1  DG1340                    1      M  44
2  DG1383                    1      F  35
3  DG1453                    1      F  55
4  DG1455                    1      F  53
5  DG1473                    1      M  33
6  DG1557                    1      M  45
7  DG1578                    1      F  57
8  DG1590                    1      M  42
9  DG1380       NON_INDUCIBILE      M  20
10 DG1381       NON_INDUCIBILE      F  52
11 DG1382       NON_INDUCIBILE      M  49
12 DG1565       NON_INDUCIBILE      M  30
13 DG1594       NON_INDUCIBILE      F  42
14 DG1620       NON_INDUCIBILE      M  40
15 DG1633       NON_INDUCIBILE      F  62
16 DG1652       NON_INDUCIBILE      M  21

```

and I filtered the Count matrix for the samples present in the design one. So I have

```
head(Counts_IGATECH_ALL_1_random_CTR)
            DG1340 DG1383 DG1453 DG1455 DG1473 DG1557 DG1578 DG1590 DG1380 DG1381 DG1382 DG1565 DG1594 DG1620 DG1633 DG1652
ENSG00000114771      0      0      0      0      0      0      0      0      0      0      0      0      0      0      0      0
ENSG00000129673      0      1      0      0      1      0      0      1      0      1      0      0      0      0      0      0
ENSG00000144452      0      0      0      0      0      0      0      0      0      0      0      0      0      0      0      0
ENSG00000143921      0      0      0      0      0      0      0      0      0      0      0      0      0      0      0      0
ENSG00000173357      0      0      0      0      0      0      0      0      0      0      0      0      0      0      0      0
ENSG00000161103      0      0      0      0      0      0      0      0      0      0      0      0      0      0      0      0

```

When I run the `DESeqDataSetFromMatrix` function I get this error:

```
converting counts to integer mode
factor levels were dropped which had no samples
Error in checkFullRank(modelMatrix) :
  the model matrix is not full rank, so the model cannot be fit as specified.
  One or more variables or interaction terms in the design formula are linear
  combinations of the others and must be removed.

  Please read the vignette section 'Model matrix not full rank':

  vignette('DESeq2')

```

and the function i ran is:

```
dds_IGATECH_1_random_CTR<- DESeqDataSetFromMatrix(countData = round(Counts_IGATECH_ALL_1_random_CTR),
                                               colData = design_IGATECH_1_random_CTR,
                                               design = ~ Age + Gender + BrS_BASELINE_PATTERN)

```

I revised my design matrix multiple times and tried to change also the combinations of the CTRs to have a different one and see if that would solve the problem but it didn't work, and I cannot figure out what it can be.

Thank you in advance for any suggestion.

Update

Design matrix:

```
dput(design_IGATECH_1_random_CTR)

structure(list(Sample = c("DG1340", "DG1383", "DG1453", "DG1455",
"DG1473", "DG1557", "DG1578", "DG1590", "DG1376", "DG1379", "DG1573",
"DG1575", "DG1579", "DG1598", "DG1627", "DG1633"), BrS_BASELINE_PATTERN = structure(c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), levels = c("1",
"2", "NON_INDUCIBILE"), class = "factor"), Gender = structure(c(2L,
1L, 1L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 1L), levels = c("F",
"M"), class = "factor"), Age = structure(c(23L, 15L, 32L, 31L,
13L, 24L, 34L, 21L, 29L, 27L, 30L, 10L, 19L, 28L, 18L, 37L), levels = c("17",
"20", "21", "22", "23", "24", "25", "26", "28", "29", "30", "32",
"33", "34", "35", "37", "38", "39", "40", "41", "42", "43", "44",
"45", "46", "48", "49", "50", "51", "52", "53", "55", "56", "57",
"58", "61", "62", "63", "64", "69"), class = "factor")), row.names = c(NA,
-16L), class = "data.frame")

```

The NON\_INDUCIBLE data are different since i have generated a new design matrix again, but the format is the same

# Answer
This typically happens when Age is a factor rather than numeric column, so each unique row is nested with itself. Design `~BrS_BASELINE_PATTERN + Gender` is fine, there is no nesting, so this narrows it down to the Age column. Cannot check as you do not provide reproducible data via `dput()`.
> 4 votes

---
Tags: r, rna-seq, differential-expression, deseq2
---


title: Algorithm for counting neighbors

text:
# Question
Title: Algorithm for counting neighbors
I need an algorithm to process protein chains. The volume size would probably not exceed 30x30x30 angstrom^3.

Say, I have a point cloud of 50,000 points in a 3D space.

I want to **count** the number of neighbors **each** point {(x1,y1,z1), (x2,y2,z2), (x3,y3,z3), ... ...,(xN,yN,zN)} has at radii {r1, r2, r3, ..., rM}.

**Which algorithm should I use?**

**Additional information:**

1. Do you need to do this count online/sequentially or all at once?
All once.
2. What restrictions might there be on the radii?
The radii will be 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0 angstrom.
3. Do you need to do this just once in this particular case or do you need a general algorithm?
I will use this to process protein sequences. Not DNA/RNA, though.
4. For a general algorithm, what's the likely range of the number of points?
The range would be 50,000 at max (although, the maximum length of available protein is 38,000 residues).
5. Are you expecting there to be relatively few or relatively many neighbors?
I am expecting relatively few neighbors, in the range of 0-50.

**NOTE:** I prefer a hash-based algorithm as it doesn't require additional data structures like KD tree, or octree.

# Answer
The usual method is dividing the space into cells, in your case 7Å cubes, assigning atoms to cells, and searching for neighbors only in 3x3x3=27 cells. This method is described in Wikipedia as cell lists, although that name is not as widely used as the method itself.

If you analyze experimental protein structures note that they may have alternative conformations (alternative locations of the same atom), and that the structure in the file represents an asymmetric unit.
> 4 votes

# Answer
K Dimensional Tree (KD-Tree) data structure might help:

https://www.biorxiv.org/content/10.1101/605907v2.full.pdf
> 0 votes

---
Tags: proteins, protein-structure, algorithms, 3d-structure
---


# Question
Title: Handling duplicates in untargeted lipidomics
I have a MS-lipidomics dataset - from human samples split into a treatment group (n = 2) and control (n = 2) - which contains several detected targets that matched to the same lipid. I don't have access to any of the raw data, only the Tier 1 and Tier 2-matched lipid names and ID numbers.

There are around 300 duplicates out of a total of 1000 detected markers, so I'm reluctant to exclude them entirely, but I couldn't find any descriptions of best-practices for what to do in this situation. I'm considering testing the correlation between duplicates and removing the ones lie below some threshold (e.g: 0.5), then averaging the remaining duplicates - but I'm not sure if this would be a valid way of approaching this.

Are there any good ways of addressing these duplicate values?

# Answer
I don't know lipids, but if 1/3 of your data set comprises duplication thats a skewed distribution with a low mean and high variance (of the duplications). It is likely to follow a negative binomial distribution. You need to plot a histogram of your data and you will find that 0 (no duplicates) is massively represented, 1 duplicate is next in frequency, 2 duplicates are lower frequency ... and on until a very small proportion of the data comprises a large number of duplicates. This is a very common distribution throughout biology.

The advice is:

* Do not remove the duplicates because they are a vital part of the distribution, i.e. biological information.
* Construct a histogram of each data set.
* Don't take "an average", i.e. don't take the mean of skewed data, thats very important if it has a high variance.

**Easy solution** The easiest way to deal with this data is take a median and construct upper and lower quartiles around the median (IQR), especially the frequency (i.e. the distribution duplicates). This will give an approximate comparison between data sets. This alone might work and certainly gives you "cut-offs".

If it is a negative binomial distribution then testing the significance between two negative binomial distributions beyond the IQR is an advanced area of statistics and a separate question.
> 0 votes

---
Tags: proteomics, multi-omics
---


# Question
Title: Gene Set Comparison Without Expression Data
I have been looking all over the web to find some answers to my problem but unfortunately, I was unsuccessful. I wish to determine whether an a priori defined set of genes in my case genes associated with Epidermolysis bullosa shows statistically significant, concordant differences between two biological states (e.g. phenotypes).

I have a list of genes from my own experiment and a predefined list of genes associated with Epidermolysis bullosa, and want to find out if the genes of Epidermolysis bullosa are overrepresented in my list.

Here is an article about GSEA, which is related to what I want to do:

https://www.pnas.org/doi/epdf/10.1073/pnas.0506580102

so I have the following data: I'm researching patients with severe or mild disease; in this research, I have found a list of 200 genes after conducting a SKAT test in R (https://www.rdocumentation.org/packages/SKAT/versions/2.2.4/topics/SKAT) with a p-value \<0.05 ( which I presume can cause a severe disease instead of mild ) I want to focus on this 200 significant genes that can perhaps affect the disease level, this genes can be -\> the L list as described in the article above. Also, the S list of predefined genes ( as is called in the article above )is associated with Epidermolysis bullosa disease. and I want to know if the genes from the Epidermolysis bullosa are overrepresented in my L list. After reading the article and many more resources, I understood that to perform GSEA (https://www.gsea-msigdb.org/gsea/index.jsp), I had to have gene expression data which i don't have .

so I wish to understand if there is a tool that could do a gene set comparison without having the expression data whether it's in R or Python that uses my gene list and predefined list ( in my case associated with Epidermolysis bullosa) ,thank you :)

# Answer
> 1 votes
You can't do a GSEA on a set alone because the differentiating statistic (or, occasionally, the rank) is used to calculate the enrichment score that is produced by the GSEA algorithm.

The other gene list type tools that I'm aware of (that my biologist colleagues use) are DAVID and PANTHER. But those are general discovery algorithms, not for comparing specific lists.

---
Tags: r, gene, statistics, gsea, go-enrichment
---


# Question
Title: Why does my HWE QQ plot have extreme deviation and what does it mean?
*This question was also asked on Reddit*

I have recently completed my thesis and one of the comments was that I report on why this graph looks this way. I have tried to find a reason but the closest I can come to is that it means that there are false positives? Or that this is possibly due to population stratification?

The trait of interest was alcohol dependency with covariates of sex and age, in a non african american population.

I am also unsure as to how this will impact my research or the implications of this?

I don't understand why this population (non-African American) "reacted" this way - do I provide a different null model? The previous populations had "relatively normal" plots. They were all run through the same pipeline.

How the data was calculated:

> There were 1061 individuals in the non-African American dataset. At QC phase 0, there were 2 individuals who had discordant sex information and 87 duplicate SNPs. When the MAF value was set to 0.000015, 125708 SNPs were removed when the MAF. Both the mind and geno parameters were set to 0.2, resulting in 127357 and 0 SNPs to be removed respectively. When setting the HWE to $1 \times 10^{-12}$, 4495 SNPs were removed. After all these parameters were set and the pipeline completed, there were 851773 SNPs left from 1044 individuals.

# Answer
> 1 votes
This much deviation from the expected distribution looks like a modeling issue, i.e. that the statistical model use to represent the null hypothesis is incorrect.

I can't see *any* values that lie along the expected line. It's more common to see a generally linear trend, with a few outliers just at the end.

It could very likely be due to population stratification. You have defined this population as "non-African American", which seems to me to be quite a vague definition that would capture a lot of different population groups with distinct genetic histories.

In my own research many years ago, I found substantially different haplotypes within the ADH gene cluster between Māori and European individuals, and I could imagine that not adjusting for population structure in other populations would have an impact on association scores.

# Answer
> 0 votes
The level of discordance is extreme, particularly for humans. This is either a publication in Nature/Science and massive advance in the field (especially humans), or the genetic factors have massive confounders.

In theory, I'm aware there must be a genetic component to alcoholism, there is a biochemical component about the rate of processing alcohol toxins. However, non-genetic (especially social/community) affects of alcoholism could predominate.

What I personally think is happening is there are very strong non-genetic confounders which are disrupting the model and this is getting confused with neutral HWE. For example, if alcoholism predominates in a community or social group the transmission is social not genetic, so the confounder is extreme because they are genetically related, but the gene(s) are not being inherited. In phylogenetic terms it will make the "gene" appear far higher up (far deeper) in the tree than it actually is. The "lone" drinker is probably the model, but group behaviour and peer-pressure is probably more realistic. Equally community effects could be suppressing alcoholism, such as religious belief.

Thus if a non-genetic model of alcoholism, (loads of social factors) resulted in a stronger fit by for example 3 logs difference, whilst the genetic model is 50 logs difference, then that is the answer and points to the need for careful population sampling for GWAS. Upvotes J Human Gen is cool.

---

As a disclaimer, I do not understand not human genetics outside immunology. I do understand multiple linear regression models and failing to account for a pivotal factor(s) causes unusual effects. An obvious example in this case is religious belief: thus alcoholism genes could be present but the environment is not presented. For example, you might conclude long hair caused a reduction in human growth because they tightly correlated in humans. However, once gender is a factored into the model, the result is clear hair-length associates with gender not height.

---
Tags: genome, snp, gwas
---


# Question
Title: How can I add a ";" at the end of each fasta header line in a file?
I have a large fasta file with RNA sequences. I need to add a `;` at the end of each header line (lines beginning with `>`) in the file.

For example, I have:

```
>Bacteria;Cyanobacteriales;Phormidiaceae;Lyngbya PCC-7419;Lyngbya sp. PCC 8106
GGAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUCUGCUUAACACAUGCAAGUCGAACGGACUCUUCGGAGUUAGU
>Bacteria;Cyanobacteriales;Nostocaceae;Nodularia PCC-9350;Nodularia spumigena CCY9414
GAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUAUGCUUAACACAUGCAAGUCGAACGGUCUCUUCGGAGAUAGUG

```

And, my desired output is:

```
>Bacteria;Cyanobacteriales;Phormidiaceae;Lyngbya PCC-7419;Lyngbya sp. PCC 8106;
GGAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUCUGCUUAACACAUGCAAGUCGAACGGACUCUUCGGAGUUAGU
>Bacteria;Cyanobacteriales;Nostocaceae;Nodularia PCC-9350;Nodularia spumigena CCY9414;
GAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUAUGCUUAACACAUGCAAGUCGAACGGUCUCUUCGGAGAUAGUG

```

# Answer
> 7 votes
One way, using and a regex address to select only the header lines and apply a substitution to append a semicolon to the end of each line:

```
sed '/^>/s/$/;/' in.fa > out.fa

```

Or, to edit the file in place:

```
sed -i '/^>/s/$/;/' in.fa

```

---

Or one way using and a ternary construct to print the line followed by a semicolon or an empty string, depending on whether or not the line is a header line:

```
awk '{ print $0 (/^>/ ? ";" : "") }' in.fa > out.fa

```

Alternatively, you may prefer:

```
awk '/^>/ { print $0 ";"; next } { print $0 }' in.fa > out.fa

```

Which is the same as:

```
awk '/^>/ { print $0 ";"; next }1' in.fa > out.fa

```

---

If you have installed, all you need is:

```
perl -pe 's/$/;/ if /^>/' in.fa > out.fa

```

# Answer
> 6 votes
... using the legendary Perl pie

```
perl -p -i -e 's/^(>.*)/$1;/' mybacteria.fa

```

**Input**

```
>Bacteria;Cyanobacteriales;Phormidiaceae;Lyngbya PCC-7419;Lyngbya sp. PCC 8106
GGAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUCUGCUUAACACAUGCAAGUCGAACGGACUCUUCGGAGUUAGU
>Bacteria;Cyanobacteriales;Nostocaceae;Nodularia PCC-9350;Nodularia spumigena CCY9414
GAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUAUGCUUAACACAUGCAAGUCGAACGGUCUCUUCGGAGAUAGUG
A

```

**Output**

```
>Bacteria;Cyanobacteriales;Phormidiaceae;Lyngbya PCC-7419;Lyngbya sp. PCC 8106;
GGAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUCUGCUUAACACAUGCAAGUCGAACGGACUCUUCGGAGUUAGU
>Bacteria;Cyanobacteriales;Nostocaceae;Nodularia PCC-9350;Nodularia spumigena CCY9414;
GAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUAUGCUUAACACAUGCAAGUCGAACGGUCUCUUCGGAGAUAGUG

```

**Details** I copy pasted your first code into a file `mybacteria.fa` then ran the Perl pie script. Perl pie is very quick, there alternative regexs (regular expressions) that would be slightly more efficient, but Perl pie is fast regardless. In fact you don't need to use regex in this example and would be faster via a `if $line == '>';` declaration (certainly true for Pythonic equivalents `if line == '>'`) but it doesn't offer the ease and flexibility, which is combined with speed, that Perl pie does. The code will work in Linux (any flavour/flavor) or Mac OSX Darwin (i.e. Terminal app)

# Answer
> 3 votes
I recommend to use Biopython as a clearer way.

```
from Bio import SeqIO

records = []
for record in SeqIO.parse("file.fas", "fasta"):
    record.description += ";"
    records.append(record)
SeqIO.write(records, "file.fas", "fasta")

```

Before:

```
>OK087550.1 Procyon lotor isolate Plo1 cytochrome b gene, partial cds; mitochondrial
TCAAACATCTCAGCATGATGAAATTTCGGCTCCCTCCTCGGAATTTGTTTGCTTCTACAGATCGCAACAG
GTTTATTCTTAGCCATGCACTATACACCAGATACAACCACAGCTTTCTCATCAGTGACCCACATTTGCCG
AGATGTAAATTATGGCTGAATTATCCGATATATACACGCTAACGGAGCTTCTATATTCTTTATATGCCTA
TTCTTACACGTAGGACGAGGCTTATACTATGGCTCCTATACATTCTCTGAAACATGAAATATTGGTGTTA
TTTTACTACTCACAGTAATAGCAACCGCTTTCATAGGTTACGTCCTCCCATGAGGACAAATATCATTCTG
AGGAGCAACAGTAATTACAAATCTCCTGTCAGCTATCCCCTATATCGGATCTAACCTTGTAGAATGAATT
TGAGGAGGGTTTTCAGTAGACAAAGCCACCCTAACACGATTCTTCGCATTCCACTTCATTCTACCATTTA
TTATCACAACGCTAGCAATAATTCACCTGCTATTCCTACACGAAACAGGATCCAATAATCCTGCTGGAAT
TACATCAGAATCTGACAAAATTCCATTTCACCCATACTACACCATTAAAGACATTCTGGGAATCCTATTC
CTTATTTTTATACTTATAGGTTTAGTGCTATTTACGCCAGACCTACTAGGAGACCCAGATAACTACACAC
CCGCCAACCCCTTAAGCACCCCACCTCACATTAAACCTGAATGATATTTTCTATTCGCCTACGCAATTCT
ACGTTCCATTCCCAACAAACTAGGAGGAGTCCTAGCCCTAGTCCTCTCCATCTTAATCCTAATCATCATC
CCACTCCTACATACCTCAAAACAACGAAGCATAATATTTCGGCCACTTAGCCAATGTTCATTCTGATTCC
TAGTAGCAGACCTCCTCGTCCTAACATGAATTGGAGGTCAACCAGTAGAATATCCCTTCATCATCATCGG
CCAACTAGCCTCCATCTTCTACTTTATAATCCTCCTGGTCCT

>OK087549.1 Procyon lotor isolate Plo17 cytochrome b gene, partial cds; mitochondrial
TCAAACATCTCAGCATGATGAAATTTCGGCTCCCTCCTCGGAATTTGTTTGCTTCTACAGATCGCAACAG
GTTTATTCTTAGCCATGCACTATACACCAGATACAACCACAGCTTTCTCATCAGTGACCCACATTTGCCG
AGATGTAAATTATGGCTGAATTATCCGATATATACACGCTAACGGAGCTTCTATATTCTTTATATGCCTA
TTCTTACACGTAGGACGAGGCTTATACTATGGCTCCTATACATTCTCTGAAACATGAAATATTGGTGTTA
TTTTACTACTCACAGTAATAGCAACCGCTTTCATAGGTTACGTCCTCCCATGAGGACAAATATCATTCTG
AGGAGCAACAGTAATTACAAATCTCCTGTCAGCTATCCCCTATATCGGATCTAACCTTGTAGAATGAATT
TGAGGAGGGTTTTCAGTAGACAAAGCCACCCTAACACGATTCTTCGCATTCCACTTCATTCTACCATTTA
TTATCACAACGCTAGCAATAATTCACCTGCTATTCCTACACGAAACAGGATCCAATAATCCTGCTGGAAT
TACATCAGAATCTGACAAAATTCCATTTCACCCATACTACACCATTAAAGACATTCTGGGAATCCTATTC
CTTATTTTTATACTTATAGGTTTAGTGCTATTTACGCCAGACCTACTAGGAGACCCAGATAACTACACAC
CCGCCAACCCCTTAAGCACCCCACCTCACATTAAACCTGAATGATATTTTCTATTCGCCTACGCAATTCT
ACGTTCCATTCCCAACAAACTAGGAGGAGTCCTAGCCCTAGTCCTCTCCATCTTAATCCTAATCATCATC
CCACTCCTACATACCTCAAAACAACGAAGCATAATATTTCGGCCACTTAGCCAATGTTCATTCTGATTCC
TAGTAGCAGACCTCCTCGTCCTAACATGAATTGGAGGTCAACCAGTAGAATATCCCTTCATCATCATCGG
CCAACTAGCCTCCATCTTCTACTTTATAATCCTCCTGGTCCT

>OK087548.1 Procyon lotor isolate Plo11 cytochrome b gene, partial cds; mitochondrial
TCAAACATCTCAGCATGATGAAATTTCGGCTCCCTCCTTGGAATTTGTTTGCTTCTACAGATCGCAACAG
GTTTATTCTTAGCCATGCACTATACACCAGATACAACCACAGCTTTCTCATCAGTGACCCACATTTGCCG
AGATGTAAATTATGGCTGAATTATCCGATATATACACGCTAACGGAGCTTCTATATTCTTTATATGCCTA
TTCTTACACGTAGGACGAGGCTTATACTATGGCTCCTATACATTCTCTGAAACATGAAATATTGGTGTTA
TTTTACTACTCACAGTAATAGCAACCGCTTTCATAGGTTACGTCCTCCCATGAGGACAAATATCATTCTG
AGGAGCAACAGTAATTACAAATCTCCTGTCAGCTATCCCCTATATCGGATCTAACCTTGTAGAATGAATT
TGAGGAGGGTTTTCAGTAGACAAAGCCACCCTAACACGATTCTTCGCATTCCACTTCATTCTACCATTTA

```

After:

```
>OK087550.1 Procyon lotor isolate Plo1 cytochrome b gene, partial cds; mitochondrial;
TCAAACATCTCAGCATGATGAAATTTCGGCTCCCTCCTCGGAATTTGTTTGCTTCTACAG
ATCGCAACAGGTTTATTCTTAGCCATGCACTATACACCAGATACAACCACAGCTTTCTCA
TCAGTGACCCACATTTGCCGAGATGTAAATTATGGCTGAATTATCCGATATATACACGCT
AACGGAGCTTCTATATTCTTTATATGCCTATTCTTACACGTAGGACGAGGCTTATACTAT
GGCTCCTATACATTCTCTGAAACATGAAATATTGGTGTTATTTTACTACTCACAGTAATA
GCAACCGCTTTCATAGGTTACGTCCTCCCATGAGGACAAATATCATTCTGAGGAGCAACA
GTAATTACAAATCTCCTGTCAGCTATCCCCTATATCGGATCTAACCTTGTAGAATGAATT
TGAGGAGGGTTTTCAGTAGACAAAGCCACCCTAACACGATTCTTCGCATTCCACTTCATT
CTACCATTTATTATCACAACGCTAGCAATAATTCACCTGCTATTCCTACACGAAACAGGA
TCCAATAATCCTGCTGGAATTACATCAGAATCTGACAAAATTCCATTTCACCCATACTAC
ACCATTAAAGACATTCTGGGAATCCTATTCCTTATTTTTATACTTATAGGTTTAGTGCTA
TTTACGCCAGACCTACTAGGAGACCCAGATAACTACACACCCGCCAACCCCTTAAGCACC
CCACCTCACATTAAACCTGAATGATATTTTCTATTCGCCTACGCAATTCTACGTTCCATT
CCCAACAAACTAGGAGGAGTCCTAGCCCTAGTCCTCTCCATCTTAATCCTAATCATCATC
CCACTCCTACATACCTCAAAACAACGAAGCATAATATTTCGGCCACTTAGCCAATGTTCA
TTCTGATTCCTAGTAGCAGACCTCCTCGTCCTAACATGAATTGGAGGTCAACCAGTAGAA
TATCCCTTCATCATCATCGGCCAACTAGCCTCCATCTTCTACTTTATAATCCTCCTGGTC
CT
>OK087549.1 Procyon lotor isolate Plo17 cytochrome b gene, partial cds; mitochondrial;
TCAAACATCTCAGCATGATGAAATTTCGGCTCCCTCCTCGGAATTTGTTTGCTTCTACAG
ATCGCAACAGGTTTATTCTTAGCCATGCACTATACACCAGATACAACCACAGCTTTCTCA
TCAGTGACCCACATTTGCCGAGATGTAAATTATGGCTGAATTATCCGATATATACACGCT
AACGGAGCTTCTATATTCTTTATATGCCTATTCTTACACGTAGGACGAGGCTTATACTAT
GGCTCCTATACATTCTCTGAAACATGAAATATTGGTGTTATTTTACTACTCACAGTAATA
GCAACCGCTTTCATAGGTTACGTCCTCCCATGAGGACAAATATCATTCTGAGGAGCAACA
GTAATTACAAATCTCCTGTCAGCTATCCCCTATATCGGATCTAACCTTGTAGAATGAATT
TGAGGAGGGTTTTCAGTAGACAAAGCCACCCTAACACGATTCTTCGCATTCCACTTCATT
CTACCATTTATTATCACAACGCTAGCAATAATTCACCTGCTATTCCTACACGAAACAGGA
TCCAATAATCCTGCTGGAATTACATCAGAATCTGACAAAATTCCATTTCACCCATACTAC
ACCATTAAAGACATTCTGGGAATCCTATTCCTTATTTTTATACTTATAGGTTTAGTGCTA
TTTACGCCAGACCTACTAGGAGACCCAGATAACTACACACCCGCCAACCCCTTAAGCACC
CCACCTCACATTAAACCTGAATGATATTTTCTATTCGCCTACGCAATTCTACGTTCCATT
CCCAACAAACTAGGAGGAGTCCTAGCCCTAGTCCTCTCCATCTTAATCCTAATCATCATC
CCACTCCTACATACCTCAAAACAACGAAGCATAATATTTCGGCCACTTAGCCAATGTTCA
TTCTGATTCCTAGTAGCAGACCTCCTCGTCCTAACATGAATTGGAGGTCAACCAGTAGAA
TATCCCTTCATCATCATCGGCCAACTAGCCTCCATCTTCTACTTTATAATCCTCCTGGTC
CT
>OK087548.1 Procyon lotor isolate Plo11 cytochrome b gene, partial cds; mitochondrial;
TCAAACATCTCAGCATGATGAAATTTCGGCTCCCTCCTTGGAATTTGTTTGCTTCTACAG
ATCGCAACAGGTTTATTCTTAGCCATGCACTATACACCAGATACAACCACAGCTTTCTCA
TCAGTGACCCACATTTGCCGAGATGTAAATTATGGCTGAATTATCCGATATATACACGCT
AACGGAGCTTCTATATTCTTTATATGCCTATTCTTACACGTAGGACGAGGCTTATACTAT
GGCTCCTATACATTCTCTGAAACATGAAATATTGGTGTTATTTTACTACTCACAGTAATA
GCAACCGCTTTCATAGGTTACGTCCTCCCATGAGGACAAATATCATTCTGAGGAGCAACA

```

# Answer
> 3 votes
```
seqkit replace -p "$" -r ";" my.fasta > my.renamed.fasta

```
* Where `-p` is the pattern to match the end of the header text `$`
* Where `-r` is the replacement pattern

https://bioinf.shenwei.me/seqkit/usage/#replace

Seqkit is an extremely powerful toolkit for working with fasta/fastq files. It can a lot with a single cli command.

---
Tags: phylogenetics, fasta, phylogeny, bash, bacteria
---


# Question
Title: DE analysis tool and method for non-coding features
I am currently working with the non-coding features of *A. thaliana*, and trying to get the DE features. Among the three DE test methods in `edgeR` such as `exactTest()`, `glmFit()` and `glmLRT()`, `glmQLFit()` and `glmQLFTest()`, which is/are the one(s) most suitable for running DE tests on miRNA, snRNA, snoRNA, pseuogenes, transposons and 5' and 3'-UTRs? Should I be using any other DE analysis tool instead of `edgeR` for some specific non-coding feature, like there's lncDIFF for lncRNA?

# Answer
> 3 votes
The post here is useful ... https://support.bioconductor.org/p/76790/

`ExactTest` is really cool and completely different from `glmQLFit` `glmFit` etc ... you'll need two controls to check it works (it could just reject everything). Its described here https://rdrr.io/bioc/edgeR/man/exactTest.html

There is a concern, some of the **lncDIFF** data looked more like a lognormal to me in the publication, but I didn't read it in sufficient detail. Anyway your data will probably approximate to a negative binomial distribution and it's statistically tricky to assess whether two negative binomial distributions are significantly different ... the `ExactTest` will do this.

Is the total output of the control significantly different to the test sample? ... the `ExactTest` will tell you. You want *p* \< 0.05, preferably a lot lower than 0.05 to be convinced. However, you need to compare two control samples against each other, where *p* \> 0.05 for the subsequent tests to be valid.

---

So really nice BUT you need to check whether it is a negative binomial FIRST. Thats what the other tests are doing. Please do check whether **lncDIFF** identified a negative binomial, because if it ain't `ExactTest` will not work, they could be using a different distribution.

The other tests you mentioned are different to the `ExactTest` because they simply assess whether the data is negative binomial. The distribution is the frequency of each RNA transcript BTW. edgeR's **lnfDIFF** requires that distribution and thats how it can assign probabilities to the difference between the frequency distribution of two identical transcripts under difference experimental conditions (or other conditions). So you need a test to say "yeah I'm a negative binomial" (assuming it is negative binomial), then everything else via **lncDIFF** is then valid.

First you've got to fit you data to a model ...

```
y <-glmGLMFit(...)

```

... is the function and that uses general linear models to fit the negative binomial. The fit will result in something called "a residual" and that is very important, the smaller the better. Testing whether that residual is small enough for a negative binomial to be declared is down to ...

```
glmLRT(y ...)

```

This is a likelihood ratio test, to test the observed likelihood against the model (p \> 0.05). You can also use

```
gof(y ...)

```

this is goodness of fit, as the link says the residual follows a chi-squared distribution, i.e. its tested on the chi2 tables and thats easy to understand. Again you need a probability \> 0.05 to declare a fit.

HOWEVER, the **lnfDIFF** guys developed a quasi-likelihood to assess the distribution

```
z<-glmQLFit(...)

```

Fits the data via quasi-likelihood, I honestly have no empirical idea how that works.

```
glmQLTest(z ...)

```

tests the fit to the distribution. You need a probability \> 0.05. Again I have not empirical idea whats going on here.

ANY non-significant result is good, doesn't matter if its GLM with gof, GLM with LRT, GLM-QL (whatever that is) because all of these are solid tests. Which would I use? The guys who developed **lncDIFF** used QL (quasi-likelihood) to assess the fit (might be a good idea). However, statistically I honestly don't know whats going on, so I would also look at the GLM and fit it against the goodness of fit - because thats easy to understand.

What do I do to assess the fit? I use machine learning, but it's way more complicated than anything here.

---

From the comments.

It's important to note the context. The OP is performing non-coding RNA work, **lncDIFF** is much better than normal **edgeR** for this type of data described here Statistical methods suitable for DE analysis of non coding RNAs.

**lncDIFF** is now available for **R** here

... so the methods of edgeR should be available to it. `glmQLFit()` is a general function in R. Providing the **lncDIFF** object will be accepted by the tests, by performing everything in **R** this should solve this.

*I would advise using* **lnfDIFF**'s QL directly, i.e.

```
y<-ZIQMLfit(...) # fit the distribution
ZIQML.LRT(y ...) # test the fit

```

The manual is here

https://github.com/qianli10000/lncDIFF/blob/master/lncDIFF-manual.pdf

---

I hope that makes sense these are pretty advanced stats so it can be tricky to understand them. Whilst I can teach code, I don't teach complex stats outside machine learning.

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


# Question
Title: Getting rid of duplicates in a dictionary
I have a tsv file that lists the reads and read lengths from a FASTA file but some reads are duplicated - that's just from the analysis I did previously - but I want to only take one instance of the read and have it be the SHORTER length than the other duplicate.

For instance, here is the beginning of the tsv file:

```
read                                    length
91290b97-5ca1-489d-b488-17afc1d217e7    3584
dad08436-3cd6-40a1-9af1-b846bd24a815    3787
f4a7bbdf-de62-48a8-a8b0-06076e251820    3637
f4a7bbdf-de62-48a8-a8b0-06076e251820    4717
21ddb2f2-654e-4124-8226-91c5a8d11a97    3786
6a2de460-a661-4c9e-ac1d-adfc3e5841e2    3520

```

You can see that `f4a7bbdf-de62-48a8-a8b0-06076e251820` is found twice with different lengths `3637` and `4717`. I want to go through and take the shorter length of that read and include it in the final tsv file, resulting in this structure:

```
read                                    length
91290b97-5ca1-489d-b488-17afc1d217e7    3584
dad08436-3cd6-40a1-9af1-b846bd24a815    3787
f4a7bbdf-de62-48a8-a8b0-06076e251820    3637
21ddb2f2-654e-4124-8226-91c5a8d11a97    3786
6a2de460-a661-4c9e-ac1d-adfc3e5841e2    3520

```

I'm trying to do this in python and have started with making the tsv file a list:

```
newlist = []
duplist = []

with open(tsv_path) as file:
    tsv_file = csv.reader(file, delimiter="\t")

    for line in tsv_file:
        if line[0] not in newlist:
            newlist.append(line[0])
        else:
            duplist.append(line)

    print("All the duplicate reads: ", duplist)
    print("The normal reads: ", newlist)

```

Which works to contain the first instance of a read in one list and the second/third instance of the read to go in a separate list with their lengths. However, the shorter length of a read is not always the first instance of a read. **Do I need to include another loop that looks at the lengths and compares each instance of a read? Is there an easier way to do this?**

# Answer
~~Not sure how you'd do this in python but~~ you need a group-by operation here. Group by the read name and pick the min length per group.

EDIT: With pandas, you can do:

```
df.groupby(['read'])['length'].agg('min')

```
> 3 votes

# Answer
You can do this with the `sort` command-line tool. This command pipe sequence sorts numerically by the second field, then chooses unique elements based on the first field

```
sort -r -k 2n,2 read_lengths.txt | sort -r -u -k 1,1 > unique_lengths.txt

```

*\[the additional `-r` bits are not strictly necessary; they help make sure that the headings stay at the top\]*
> 2 votes

---
Tags: python, file-formats, pandas, duplications
---


# Question
Title: How to rearrange gene ontology obo file?
I have a gene ontology obo file, and would like to rearrange the obo format file in a single column to a text file in three columns comprised of 'parent (id:GO:xxxxxxx)', 'child (is\_a: GO:xxxxxxx, part\_of GO:xxxxxxx, or positively/negatively\_regulates GO:xxxxxxx. For example...

Thank you for your help in advance!

Original file:

```
[Term]
id: GO:0000001 # parent
name: mitochondrion inheritance
namespace: biological_process
def: "The distribution of mitochondria, including the mitochondrial genome, into daughter cells after mitosis or meiosis, mediated by interactions between mitochondria and the cytoskeleton." [GOC:mcc, PMID:10873824, PMID:11389764]
synonym: "mitochondrial inheritance" EXACT []
is_a: GO:0048308 ! organelle inheritance # child
is_a: GO:0048311 ! mitochondrion distribution # child

[Term]
id: GO:0000015 # parent
name: phosphopyruvate hydratase complex
namespace: cellular_component
def: "A multimeric enzyme complex, usually a dimer or an octamer, that catalyzes the conversion of 2-phospho-D-glycerate to phosphoenolpyruvate and water." [GOC:jl, ISBN:0198506732]
subset: goslim_metagenomics
synonym: "enolase complex" EXACT []
is_a: GO:1902494 ! catalytic complex # child
relationship: part_of GO:0005829 ! cytosol # child

[Term]
id: GO:0000122 # parent
name: negative regulation of transcription by RNA polymerase II
namespace: biological_process
alt_id: GO:0010553
alt_id: GO:0045816
def: "Any process that stops, prevents, or reduces the frequency, rate or extent of transcription mediated by RNA polymerase II." [GOC:go_curators, GOC:txnOH]
synonym: "down regulation of global transcription from RNA polymerase II promoter" RELATED []
synonym: "down regulation of transcription from RNA polymerase II promoter" EXACT []
synonym: "down-regulation of global transcription from RNA polymerase II promoter" RELATED []
synonym: "down-regulation of transcription from RNA polymerase II promoter" EXACT []
synonym: "downregulation of global transcription from RNA polymerase II promoter" RELATED []
synonym: "downregulation of transcription from RNA polymerase II promoter" EXACT []
synonym: "inhibition of global transcription from RNA polymerase II promoter" RELATED []
synonym: "inhibition of transcription from RNA polymerase II promoter" EXACT []
synonym: "negative regulation of gene-specific transcription from RNA polymerase II promoter" RELATED []
synonym: "negative regulation of global transcription from Pol II promoter" RELATED []
synonym: "negative regulation of transcription from Pol II promoter" EXACT []
synonym: "negative regulation of transcription from RNA polymerase II promoter" EXACT []
synonym: "negative regulation of transcription from RNA polymerase II promoter, global" RELATED []
is_a: GO:0006357 ! regulation of transcription by RNA polymerase II # child
is_a: GO:0045892 ! negative regulation of DNA-templated transcription # child
relationship: negatively_regulates GO:0006366 ! transcription by RNA polymerase II # child

```

Desired format:

```
Parent       Child       Relationship
GO:0000001   GO:0048308  is_a
GO:0000001   GO:0048311  is_a
GO:0000015   GO:1902494  is_a
GO:0000015   GO:0005829  part_of
GO:0000122   GO:0006357  is_a
GO:0000122   GO:0045892  is_a
GO:0000122   GO:0006366  negatively_regulates

```

# Answer
One potential solution using awk:

```
awk 'BEGIN {
    OFS = "\t"
    print "Parent", "Child", "Relationship"
}

/GO:[[:digit:]]{1,}/ {
    if ($1 ~ /^id:/) {
    parent_id = $2
} else if ($1 ~ /^is_a:/) {
        gsub(":", "", $1)
    print parent_id, $2, $1
} else if ($1 ~ /^relationship:/) {
        print parent_id, $3, $2
    }
}' test.txt
Parent      Child       Relationship
GO:0000001  GO:0048308  is_a
GO:0000001  GO:0048311  is_a
GO:0000015  GO:1902494  is_a
GO:0000015  GO:0005829  part_of
GO:0000122  GO:0006357  is_a
GO:0000122  GO:0045892  is_a
GO:0000122  GO:0006366  negatively_regulates

```

---

Including comments:

```
awk 'BEGIN {                      # before processing the first line
    OFS = "\t"                    # set the output field separator to tab
    print "Parent", "Child", "Relationship"    # and print the header
}

/GO:[[:digit:]]{1,}/ {    # for any line that contains a GO term
    if ($1 ~ /^id:/) {            # if the line starts with id:
    parent_id = $1                # set the parent id
} else if ($1 ~ /^is_a/) {        # then, for lines that start with is_a
        gsub(":", "", $1)         # remove the colon
    print parent_id, $2, $1       # and print the parent, child and relationship
} else if ($1 ~ /^relationship:/) {  # then, for lines that start with relationship
        print parent_id, $3, $2   # print the parent, child, and relationship
    }
}' test.txt

```
> 2 votes

---
Tags: gene, gene-ontology
---