# Question
Title: How do I color selective residues on the surface of a protein in pymol?
I am somewhat familiar with the pymol userface, but not very good with commands. I'm working with a pentameric subunit protein and I want to color specific residues on the surface of the protein, not just the cartoon. I can use select resi #-#, and color those with the (sele) tab on the right, but I do not know how to color the surface this way. When you select the sequence it will only color the cartoon. Any suggestions?

I'm using PyMOL 2.3.5 on Windows.

# Answer
> 1 votes
To clarify, you are seeing reds and blues on the surface? If so, these are the sidechain oxygens and nitrogens. So you need to colour the whole residue not only the `element C` atoms. `color 0xFF7F50, resi 123` will colour the residue in the cartoon and surface as I omitted and `element C`. To revert `color atomic`.

In PyMOL an atom has a colour, which is then used by the representations. To have a differentially coloured representation, you need to make a copy of the object. `create copied, obj_name` and `hide cartoon` and `show surface` will do it.

---
Tags: pymol
---


title: Find SNPs in yeast genomes

text:
# Question
Title: Find SNPs in yeast genomes
I'm a new Bioinformatic scientist working for a yeast genetics company.

**Objective** To create a database of yeast genomes from NCBI and identify SNP variants.

In my pipeline

1. FastQC,
2. Trimmomatic,
3. BWA
4. GATK

The method being to check the quality of the downloaded reads, trim off any unnecessary information, and to align them prior to VCF.

NCBI sequence data comes from both *Illumina* sequencing and *PacBio*.

**Questions**

1. Is the pipeline reasonable, for example, is BWA ok over more recent developments?
2. Could the same pipeline be used both for Illumina and PacBio sequences?

**Note** I don't plan on assembling any of the reads I have as I'm just looking to comb through to find variations.

# Answer
> for illumina reads (only)
>
> prerequisite: download n prepare

> the fna n gff file of yeast reference

> mapping\_in\_02.directory

1st step is mapping the downloaded reads to the reference genome:

1. build index for the reference using the following command

bwa index ../01./GCF\_000146045.2\_R64\_genomic.fna

> index the reference genome for rapid searching and aligning

2. do the mapping using bwa mem (better with bash shell script, then there will be no need for each individual one by one)
3. After mapping, you could use GATK4 for the variation calling based on the mapping result.

The steps of variation calling by GATK4 include dealing with the mapping result ( mark duplicates and add read group ), raw variation identification ( joint call followed by single individual variation calling ), quality recalibration and realign, second round of variation calling and final variation filtering ( filtering is applied to haplotyping to remove uninformative mutations. An obvious filter is to remove SNPs with identical calls in all the samples ).

The final variation calling result with SNPs and InDels could be obtained and the final vcf file could be sent to SnpEff for further annotation ( an easier way is using SnpEff in Galaxy with the right reference genome version, eg R64-1-1.75 ).

The isolates from NCBI would indicate that there are 16 chromosomes and one mitochondrial chromosome after SnpEff analysis. In the analysis, the numbers of insertions and deletions will be obtained.
> 3 votes

---
Tags: snp, variant-calling, genetics
---


title: change/edit FASTA headers

text:
# Question
Title: change/edit FASTA headers
I have a FASTA file with 97,000 entires arranged like this :

```
>tr|A0A075B6G3|A0A075B6G3_HUMAN Dystrophin OS=Homo sapiens OX=9606 GN=DMD PE=1 SV=1

```

I want to have an output like this:

```
>A0A075B6G3|1|training

```

I am aware of many one liners in `sed`, `perl` and `awk` but none are working . I have something like this :

```
awk '/>/{$0=">ID|1|training"}1'

```

Kindly help.

# Answer
> 7 votes
One way using :

```
awk 'BEGIN { FS=OFS="|" } /^>/ { print ">" $2, "1", "training"; next }1' file.fa

```

This sets the *field separator* (FS) and the *output field separator* (OFS) to a pipe character. Then for each header line, we can just print the fields we want. The output fields will be delimited with pipe characters since we already set `OFS` to a pipe character in the `BEGIN` block. We can use *next* here to force awk to immediately stop processing the current record and move on to the next record, which is just the next line, since, by default, awk's *record separator* is a newline character.

Since the `1` on the end just returns true, it performs the default action, which is to print the current record stored in `$0`. It is just the same as:

```
awk 'BEGIN { FS=OFS="|" } /^>/ { print ">" $2, "1", "training"; next } { print $0 }'

```

Which is the same as:

```
awk 'BEGIN { FS=OFS="|" } /^>/ { print ">" $2, "1", "training"; next } { print }'

```

# Answer
> 1 votes
```
perl -p -i -e 's/>[a-z]+?\|([A-Z]+[^\|]+?\|).*SV=([0-9]+)/>$1$2\|training/g' myfilename.fa

```

Not tested yet. Just check the sequence ID starts with a capital A-Z character. Probably wise to test it without the `-i`

# Answer
> 1 votes
Here's another simple perl solution:

```
$ perl -lpe 's/\|.*(\d+)/|$1|training/ if s/>.+?\|/>/' file.fa
>A0A075B6G3|1|training

```

Note that you didn't tell us where the `1` came from, so I am taking the last numbers on the line, the `SV=1` in your case. If you want to hard-code the `1`, you can do so easily with:

```
perl -lpe 's/\|.*/|1|training/ if s/>.+?\|/>/' file.fa

```

The `-lpe` options mean "**p**rint each line, adding a new**l**ine character to each `print` call, after applying the script given by `-e`". So this will read the file line by line, do its thing and print the result.

The `s/old/new/` is the substitution operator and will replace `old` with `new`. In this case, we are replacing a `|` character (`\|`) and everything after it (`.*`) until the last set of numbers on the line which we capture and save as `$1` with `(\d+)`. This substitution will only happen if we have already successfully substituted the leading `>` and everything until the first `|` (`>.+?\|`) with just a `>`.

If you want to edit the original file in place, use the `-i` flag (but give it a value, `.bak` in the example below, since that ensures a backup file with the original data and that extension will be created):

```
perl -i .bak -lpe 's/\|.*(\d+)/|$1|training/ if s/>.+?\|/>/' file.fa

```

---
Tags: fasta, awk
---


# Question
Title: Is there an alternative to bulked segregant analysis for insects?
This strategy seems to be most commonly used for plants. When crossing animals isn't possible, how can I do a similar study to identify a particular locus responsible for a polymorphic trait (specifically, the pigmentation of a beetle)?

# Answer
I'll have a go at answering this. I don't understand beetles within Insecta: I am well aware they very distinct, obviously in phenotype and ecology, but also genetically. I think they have unique compositional bias even amongst holometabolous insects.

Personally I think if you've a recent genetic trait, e.g. insecticide resistance, that evolved within the last 30 years or less, you're dealing with pretty homogeneous backgrounds between traits. In this example trait, bulk segregant analysis might work.

Beetles are really ancient. I suspect their families are particularly ancient, because they are just so successful and I'd guess Scarabaeidae are ancient too, i.e. the family scarab beetles. This family could easily be 50 million years old, or older. Each species therein could a few million years old. For insects thats normal and these are beetles.

The concern here is over million+ years of evolution (I don't know the study size), they will comprise loads of neutral signal that would act as noise, i.e. there will be excessive third codon wobble which will mask SNPs/genes causing colour (color) variation.

My approach, given I know nothing about colour (color) variation in insects other than "white-eyed" phenotype in Drosophila melanogaster,

1. assess genetic complexity of colour variation in other insect Orders;
2. get hold of as many proteins involved in colour in insects as possible;
3. particularly get hold of colour variation proteins in beetles;
4. sequence the genome of a representative scarab beetle (which I think your doing)
5. perform `psi-blast` using the hypotheses from 1., 2. and particularly 3. to assess which proteins in the scarab beetle might be responsible.
6. perform deep sequencing on those targets, because they could be multigene families *****.

The next step is phylogenetics to assess convergent and/or parallel evolution **\****. That in my view is a separate question because its quite complicated and you'd need to write your own pipeline to automate it (you could do it by hand however).

What I would be cautious about is that beetles are not any other insect so inter-Order homology needs to be treated with caution.

---

***,** Just to point out, there are lots sequencing strategies its not my thing, e.g. low coverage draft genomes. PCR deep sequencing for a given gene target, without any draft genomes and isn't conserved e.g. rRNA types genes, will need small PCR fragments because of the neutral mutations which will disrupt primer specificity. If you're good at sequencing a draft genome could be cool (just guessing).

\****,** A phylogenetics approach would dovetail with the COI and rRNA gene sequencing (28S?) already underway, because these genes should approximate to the species phylogeny. What is sought are residues within the colour genes (or perhaps presence/absence genes) that will cluster not according to their evolutionary descent but according to their colours (colors).
> 2 votes

---
Tags: wgs, genotyping
---


# Question
Title: How does one distinguish nuclear DNA from mitochondrial DNA when doing WGS?
I'm interested in doing de-novo sequencing but also phylogenetic analysis. In particular, after de-novo sequencing and annotating the genome, I need to align the CO1 gene and the nuclear 28S rRNA gene of several species. When extracting DNA and sequencing with PacBio for example, how does the assembler know what corresponds to nuclear and what to mitochondrial DNA?

# Answer
According to this pdf, from pacb.com

isolation with mitoISO MITOISO

The kit will include three steps: 1) mechanical rupturing of tissues or cells 2) removing cellular debris and nuclei by low speed centrifugation and 3) harvesting mitochondria by high speed centrifugation.

Then, you can follow the DNA extraction with mitochondrial DNA, and then sequenced on the PacBio Sequel.

If you are doing whole genome sequencing, the mitogenome will have a higher read depth according to this article stating that lower read depth of the nuclear genome separating the two different origins:

> those of mitochondrial origin & those of nuclear origin
> 1 votes

---
Tags: assembly, phylogenetics, mtdna
---


# Question
Title: Can we do NGS library preparation using UMI with large amount of DNA input?
We know that in next-generation sequencing (NGS), the unique molecular identifier (UMI) can reduce or eliminate sequencing or PCR errors and result in very high accurate data. Therefore, UMI is widely used in rare mutation detection. A typical application is cancer early detection or minimal residual disease (MRD) detection that try to find the very rare DNA from cancer cells in the blood.

However, we also know that the cell-free DNA in the blood is very low, typically 5 ng in 1 mL blood. And it usually takes about 10 mL blood to get about 50 ng DNA to get sequenced with UMI. Limited DNA input leads to many false negatives.

So I came up with an idea: **What if I got a sample with plenty of DNA? Can I detect rare mutations using UMI in this sample?**

The UMI library preparation protocols I have met all use a small amount of DNA input, say 50ng. **If I have 1μg DNA, do I have to cut it to 50ng to prepare the UMI library?**

# Answer
The question is

> "Do I add 20-fold more DNA than the kit states?"

Its not my area - *this is wet-lab* \- but there will be a tolerance however 20-fold is excessive (\>1 log), i.e. no it doesn't help and could cause (expensive) downstream problems. The people to contact are the support for the given kit, there are unlikely to support your idea.

The answer isn't wet-lab but almost certainly is bioinformatics and detecting low frequency somatic driver mutations. The Broad Institute worked on this problem for years, its nothing new, i.e. identifying driver mutations in highly mixed sample comprising healthy somatic tissue and cancer cells.

I've tracked their work and have a reasonable understanding of what they do and why: oncology isn't my field. Your in-house bioinformatician should know this stuff.

One suggest is using MuTect from the Broad Institute, I assume there have been many custom solutions since then. Their famous product is GATK for VCF, I don't think thats the solution here.

Again cancer isn't my area BTW this is just common knowledge within bioinformatics, i.e. knowing what someone else is doing and why is part of the knowledge base.

Please do try to retain a bioinformatics focus for future questions.
> 1 votes

# Answer
Whole genome sequencing uses between 50 ng and 1 μg of high-quality DNA (A260/A280: 1.7–2.0) for PCR-amplified or non-amplified libraries if your focus is precision oncology as described in the following open access review in Seminars in Cancer Biology paper by Meggendorfer *et al* 2022
> 1 votes

# Answer
Do you mean the following type of technology?
> -1 votes

---
Tags: ngs, mutations, library, umi
---


# Question
Title: Fixing FASTA file for Local BLAST Database
I recently prefetched 157 SRA files from an NCBI BioProject using the SRAtoolkit. I then used the toolkit to download those files in FASTA format.

Each individual FASTA file looks something like this:

*(example for SRR5678544 file)*

```
>SRR5678544.1 HWI-ST1146_88:5:1112:6472:81473 length=194
CCATGCGGGGTATCGTATGCTTCCTTCTGCACTACCCTTTAGCTGTTCTATATGCTGCCACTCCTCAATTGGATTAGTCTCATCCTTCAATGCTATCACAAGAAGTAGAGACACAAATGCAAGAGGAGCATATAAATTACAAAACACCATCACTGAGGGCCCTAAAGCGGTTCCCACGAAAAAAAGGAGAGTAG
>SRR5678544.2 HWI-ST1146_88:5:1113:13218:62635 length=194
 TGGATGGTGGTTGGGAGTGGTAAGGTTGAATGAGACACGGTAACGAGTGGGGAGGTAGGGTAATGGAGGGTAAGTTGAGAGACAGGTTCGTCAGGGGACACACCACACCACACCCACACCCACACACCCACACACCACACCCACACACCACACCCACACCCACACCCCCACACACCACACCCCCACCCCACCCA
>SRR5678544.3 HWI-ST1146_88:5:2206:11224:22269 length=194
 GTCTCTTAACTTACCCTCCATTACCCTACCTCCCGACTCGTTACCCTGTCTCATTCAACCATACCACTCCCGACCACCATCCATCCCTCTACTTACTACAGTATGGTGAGTGGGACATGGTGGATGGTAGGGTAAGTGGCAGTGGAGTTGGATATGGGCAATTGGAGGGTAACGGTTATGGTGGACGGGGGGTG

```

And so on and so forth. (Ignore the " and ' marks I had to input those to format it right on here)

**Question 1:** Can someone explain why this file has so many sequences (\>)? Every FASTA file I have worked with has one line with "\>" at the top and a header sentence followed by a single big paragraph of the genome sequence.

**Question 2:** If I wanted to combine all the sequences in a single FASTA file to look like a "normal" FASTA file - just a single paragraph - how would I do that?

**Reasoning:** I'm asking simply because I need to make a BLAST database out of the 157 sequences and the "makeblastdb" command only takes in a single file as the input. I was going to try and combine each individual FASTA file into single paragraphs and then make a huge FASTA file that includes all 157 sequences separated by "\>".

# Answer
> 2 votes
Briefly the best strategy I think is to download fastQ files and process this via an assembler. I wouldn't enter raw short read Illumina data into a local Blast database. What I'm pretty certain of is thats a bad move.

Which assembler I personally think is a separate question on Bioinformatics SE. It will depend on the genome, for large eukaryote genomes BWA-mem2 seems popular, but its not my thing and there's other stuff upstream and downstream (someone might dispute that). For assembling small genomes, Spades is popular, but preferences are changing and there's also upstream QC. The other question is whether you have a reference genome to assist assembly.

---

The update is this a yeast genome with a reference. This is complicated because its a small genome but will comprise heterozygosity. If it was me - i.e. someone way outside the yeast community - I'd use Spades and provide a reference genome to perform an assembly. However, I'm not sure thats at all cool, because you will will reduce the eukaryotic genome to a haplotype. This could get complicated because what you would be left with may not be biologically real, e.g. two alleles at different parts of the same chromosome on different strands may result in a *in silico* artefact. I think its a separate question personally.

---
Tags: phylogenetics, fasta, phylogeny, makeblastdb, sra
---


# Question
Title: correlation between imputed genotype and true genotype
There is a graph in the link https://odelaneau.github.io/GLIMPSE/rsquare\_eur.html. This graph shows the imputation accuracy of certain MAF bins. Now I understand that $r^2$ is the correlation between imputed genotype and true genotype. What I don't understand is how they calculate this $r^2$ correlation. Can anyone suggest to me how I can define the r^2 correlation between true genotype and imputed genotype? I suppose it depends on the individual who defines his own $r^2$. So, it may differ for the different algorithms. However, different aspects or different $r^2$ definitions will help me understand and represent this graph slightly differently.

# Answer
This is done by downsampling. Take the 1000 genomes, set some genotypes as missing `./.`, impute them using GLIMPSE, then measure correlation between the genotype dosages in the imputed dataset and the same sample at high coverage - in this case, the 1000 genomes sequenced to 30x coverage by the NY Genome Centre.
> 2 votes

---
Tags: correlation, imputation
---


# Question
Title: What is the best way to process yeast genomes?
I have obtained several hundred raw, unassembled yeast genomes from NCBI and I am looking for advice on how to process the genomes for downstream analysis.

I have a reference genome (S288C) to use for assembly/alignment but I'm not sure what to do first.

So far I have run a quality assessment (FastQC) followed by trimming (Trimmomatic). But I'm not sure if I should align or assemble next. The problem with yeast is if I assemble to a reference genome it reduces the genome to a haplotype and could miss moments of heterozygosity present. And I'm not sure how it would work to do a de novo assembly as I'm assuming it's the same problem. Same with aligning to the reference genome.

Are there ways to work around this? Does the yeast research community just ignore this issue? Are there tools currently out there that take all of this into account?

# Answer
> 2 votes
What is your data type? Many assemblers can resolve heterozygosity.

PacBio HiFi reads can be used with hifiasm to resolve both haplotypes rather trivially. This is the best case scenario and older sequencing projects are unlikely to have HiFi reads.

If they are merely PacBio subreads, you can use Falcon-Unzip to get heterozygous regions (though this is now quite old and I am not sure that it is recommended anymore).

For Illumina, the options are worse. Platanus is an assembler intended for the purpose of assembling heterozygous genomes. If nothing else, you can make the assembly and then just call variants.

**Note that many yeast are probably sequenced from a haploid strain.**

---
Tags: sequence-alignment, genome, assembly, genotyping, haplotypes
---


# Question
Title: Single or partitioned model for supermatrix tree inference?
I am dealing with several genes in my dataset. For each of the genes, I have built gene tree, with their best estimated model. I am intended to look into the effect of the tree from concatenated sequences (supermatrix tree), using Maximum Likelihood approach.

I am considering two options to apply the model for the supermatrix tree inference.

1. Obtain a single best estimated model for the supermatrix gene
2. Apply partitioned model for the supermatrix gene, based on each gene's model

From what I have read in most if not all of the research articles, they seemed to be using option 1, where a single model applied for the supermatrix tree inference. Since the purpose of using evolutionary model is to describe the changes in the gene and given that the evolution pattern in each gene differ, why wouldn't the respective model of each gene be considered in the tree inference?

For examples,

Dataset A: gene1 and gene2 are having LG and WAG model respectively as the best estimated model. But when the two genes are being concatenated, it turned out that LG is the best model for the matrix. It seemed that most if not all of the experiment opts for option 1 in this case. Theoretically, will the supermatrix tree using option 1 and option 2 in the model application still similar? Also, will the evolutionary pattern in gene2 still being represented correctly on the tree when we opt for option 1 in the supermatrix? Or it actually does not impact much on the calculation of the tree?

Dataset B: gene1 and gene2 are having LG as their best estimated model. The concatenated matrix has LG as its best estimated model too. Will model partitioning has any effect in this case? Does the ML tree/site likelihood calculation affected by the length of the input sequences?

# Answer
Okay, you can't concatenate genes into a singular matrix particularly if these are rRNA genes versus protein genes, except in rare situations where the tests say 'thats okay' (below). This is simply because they evolve at very different rates rRNA genes have loci which evolve very slowly and are highly conserved.

The theory goes back to incongruence length difference test (ILD), which is somewhat confusingly also known as the partition homogeneity test and was then modified to become the localised ILD test.

The modern version of this test, simply called 'homogeneity' is found in IQ-Tree. Thus, the best approach is to test the assumption empirically,

> IQ-TREE provides three matched-pairs tests of symmetry (Naser-Khdour et al., 2019) to test the two assumptions of stationarity and homogeneity. A simple analysis:
>
> iqtree2 -s example.phy -p example.nex --symtest-only

Just to be clear this tests whether the rate between two partitions (genes) is homogeneous. If 'yes' the supermatrix is fine, if 'no' the data needs partitioning.

The theory is simply that if two genes evolve at noticeably different rates they can't be combined into a singular matrix. The theory goes back a long way is used to circumvent long-branch attraction issues. Thus, at a minimum the rRNA genes and protein genes need to be compartmentalised, i.e. have their own matrix. Between protein genes, if they are house-keeping genes, generally they'll have homogeneous rates (which you can empirically test via IQ-TREE).

If someone is concatenating rRNA with protein genes without accommodating rate variation via partitioning it is simply sloppy phylogenetics and if published represents a shortcoming with the reviewing process. This is because in all certainly these loci will fail any rate homogeneity test, there's a slim off-chance its okay which you can test it empirically. Say if 18S and 28S were concatenated (or just 16S if these are bacteria) with an atp gene, that can give a different tree to 18S with atpA, atpB, atpC, atpD etc... so the ratio between rRNA to protein genes can start influencing the tree ... which is simply wrong.

If someone concatenates a load of house-keeping genes thats okay because there's a good chance they're homogeneous and if even the test failed, they may not be that far out, rRNA versus proteins can be way out. This can easily be an order of magnitude out and thats difficult to ignore. When there's proteins under selection versus conserved housekeeping genes ... its very unlikely homogeneity would hold.

**Summary** As a general rule of thumb in your method of tree building compartmentalising all genes is safer route, computationally there's no difference between zero partitions and a partition for every gene. You are never going to make a mistake.

Where it gets complicated is when the phylogenetic parameters, particularly in advanced models, are themselves estimated via maximum likelihood alongside the topology. The rational for this is complicated. Computationally, this can get very, very expensive even on supercomputers and hence there's a tendency to make an assumption (often without checking it) that there are fewer partitions than there actually are. Regardless, assuming away rRNA versus protein genes even in this method of tree estimation should not done and the test will almost certainly fail.
> 0 votes

---
Tags: phylogenetics, phylogeny, modelling, substitution-model
---


# Question
Title: Where can I download 30x 1000 genomes cram files?
From the preprint published by 1000 genome project (https://www.biorxiv.org/content/10.1101/2021.02.06.430068v1.full) I think the 30x data is for WGS. Can anyone confirm for me if the following file is really WGS data? ftp://ftp.sra.ebi.ac.uk/vol1/run/ERR324/ERR3242643/HG03007.final.cram.
Also, Is there any more freely available sequencing data where I can find WGS data for 30x and low coverage(\<10x) for the same sample for the South Asian population?

# Answer
Yes, those are the .cram files for the 30x whole genome sequencing data for 1000 genomes.

If you want 10x and 30x sequence data for the same individual, you can just downsample the 30x .cram file using `samtools view -s`.
> 2 votes

---
Tags: wgs, 1000genomes
---


# Question
Title: Methylation data: beta values are normally distributed
I am trying to incorporate Methylation data (Illumina 450K) into my pipeline. I was provided with two versions of the same data: one (dataset $a$) normalized with ComBat and the other (dataset $b$) normalized using the residuals method (i.e., taking the residuals of the models fitted on the covariates). Looking at the distribution of the beta values, I noticed that:

* The values for dataset $a$ have a typical bimodal distribution with peaks at around $0$ and $1$
* The values for dataset $b$ are roughly normally distributed with the mean at around $0.5$.

I am not that familiar with Methylation data, so maybe this is expected. Is it? I did not find any paper showing similar results.

# Answer
$\beta$ values are the ratio of methylated to unmethylated probe intensities. This means that $\beta \in \[0, 1\]$, so it can't be normally distributed, by definition. Its distribution may have sections that look like bell curves, but it isn't normally distributed. You can take a look at this paper by Du *et al.* \[1\] for details about the statistics of $\beta$ values.

That said, the distribution of all $\beta$ values from an array or bisulfite sequencing (or some other DNA methylation assay) should be bimodal with modes at 0 and 1. Nearly all CpGs in the human genome are methylated or unmethylated. There may be some hemi-methylated sites, but for the 450K array you mentioned, there should be few of these.

It sounds like you might be looking at some residual value, or some other value from the analysis in the paper for dataset $b$ other than $\beta$. Could you be looking at $M$ values, or something else?

\[1\]: Du *et al.*, BMC Bioinformatics, 2010
> 2 votes

# Answer
It's not clear from your question what the statistical design of the experiment is, so I can't judge whether or not a distribution around 0.5 would be expected from the residuals after fitting.

Given that dataset *a* looks like what you expect, I'd be inclined to use that (and perhaps ask the person who gave you this data why dataset *b* looks so different, and why it ended up with a 0.5 mean).
> 0 votes

---
Tags: methylation, multi-omics, combat
---


# Question
Title: Is loss/gain of function reflected in RNA-seq transcript counts?
Do LoF/GoF transcripts count toward the RNA-seq TPM count? Or would these LoF/GoF transcripts only be detected by isoform quantification?

# Answer
# Short answer

Since the effects of any given mutation is so context-dependent (and deciding which effects matter will change between researchers), the process of "variant annotation" is typically separate from RNA-seq mapping/quantification.

If you use a common RNA-seq mapping/quantification tool like STAR, HISAT2, or Kallisto, these counts will combine transcripts with and without the mutation of interest. This is because for the purpose of mapping, mutations are no different than sequencing errors, and when you're mapping you don't actually know what is what.

# Long answer

You might not be able to quantify GoF/LoF at the isoform level, either. Again, this is because the effects of mutations are so context-dependent.

If you want to know how many transcripts do or don't contain a specific mutation, you'll have to work backwards from that mutation and look at the RNA-seq reads themselves. You can use tools like Rsamtools, htseq-count, or even IGV for this.

If you want to know whether a given protein's function is affected, you have to consider all the types of mutations than can lead to that protein being affected (in your biological scenario of interest, of course).

* Some proteins may lose their binding affinity if a single amino acid is changed, in which case you can look for non-synonymous mutations in the RNA-seq reads.
* Some genes may undergo alternative splicing and produce different proteins because of some mutations. You could find this with a differential isoform expression analysis.
* Some genes could be deleted on one allele, leading to a lower expression of that protein. Biologically, that protein is still expressed, but having less of it makes it less functional in general. You might be able to find that with copy number estimation or allele-specific expression analysis.

If you're looking at non-coding RNA, what mutations affect their function will be different from what you'd look for if you care about protein function.

All in all, any investigation about what a GoF/LoF mutation looks like and how to go about finding it will depend on you biological scenario. And you may not find it by just looking at the TPM counts of a gene or its isoforms.
> 1 votes

---
Tags: rna-seq, gene-expression, rna, rna-alignment, rna-structure
---


# Question
Title: How to quantifiy of specific genes from shotgun metagenome?
I have googled a "lot", couldn't find any specific answer to the question. So, I am here seeking for your guidance. My question is similar to this. I have several metagenome (n=30). But for simplicity’s sake, let's assume only one metagenome. I wanted to quantify lasA genes in the metagenome. So far, what I have done (not sure if it's the right approach) is:

1. Download the lasA genes from uniport.
2. I then made the multiple sequence alignment (lasA.msa) of the lasA gene downloaded in the previous step.
3. The protein profile was generated `hmmbuild res.hmm lasA.msa`
4. Genome assembly (contigs) were done with megahit `megahit -1 r1.fastq.gz -2 r2.fastq.gz --min-contig-len 1000 -m 0.85 -o ASSEMBLY/ -t 20`
5. ORFs were predicted and extracted using `getorf` from the contigs (from step 4)
6. The profile (from step 3) was used to search the ORFs using `hmmsearch`
7. Several hits were obtained from `hmmsearch` indicating that the lasA gene is present in the metagenome

Now, my question is, how do I count the number of reads affiliated to the lasA gene? I feel that I need to map the reads to the assembly, but I don't know how to go about that.

NOTE: I have used MG-RAST, Humann3, and eggNOG-mapper, but none of them showed the presence of lasA gene, however the above-mentioned hmmsearch was able to find the gene, what could be the reason? Problem with my approach? Or hmmsearch using the protein profile is more sensitive?

Thanks and regards

# Answer
> 1 votes
The way I would do it is to make an assumption that the reads mapping to your contig == (abundance of the genes in the contig)/contig length.

1. Map raw reads to the contig containing your hit (original contig, pre-translation)
2. Normalise it with RPKM:

RPKM = NumberOfReadsThatAligned / ( ContigLength/1000 * totalNumberReads/1,000,000 )

3. This will give you the relative abundance of the contig which is an APPROXIMATION of the abundance of your gene.

Alternatively, to be more precise you can FragGeneScan your raw reads, extract the sequence of your putative protein and map translated raw reads to it.

For your other question: Yes, HMMs are more sensitive than blast or diamond since its output is probability score that these protein contain motifs extracted from your alignments (given that there is one). For that reason you need to look at the bitscore and E-value. If the bitscore is low and/or E-value is high, you have a high chance of false-positive.

Best wishes,

A

---
Tags: genome, assembly, metagenome, bowtie2, hmmer
---


# Question
Title: MarkDuplicatesSpark failing with cryptic error message. MarkDuplicates succeeds
I have been trying to follow the GATK Best Practice Workflow for 'Data pre-processing for variant discovery' (https://gatk.broadinstitute.org/hc/en-us/articles/360035535912).

This has all been run on Windows Subsystem for Linux 2 on the Bash shell.

I started off with FASTQ files from IGSR (https://www.internationalgenome.org/data-portal) and performed alignment with Bowtie2 (instead of BWA).

This produced a SAM file which I converted to a BAM file using SAMTOOLS.

I have been trying to use MarkDuplicatesSpark (https://gatk.broadinstitute.org/hc/en-us/articles/4409897162139-MarkDuplicatesSpark) but am getting an output that I can't make sense of.

I am running the following command:

```
gatk MarkDuplicatesSpark -I HG00102_hGRCH38_exome_aignment.bam -O HG00102_hGRCH38_exome_aignment.marked_duplicates.bam --remove-sequencing-duplicates --spark-master local[*]

```

I am getting the following output:

```
Using GATK jar /mnt/c/Users/angus/Documents/Bioinformatics/GATK4/gatk-4.2.4.0/gatk-package-4.2.4.0-local.jar
Running:
    java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /mnt/c/Users/angus/Documents/Bioinformatics/GATK4/gatk-4.2.4.0/gatk-package-4.2.4.0-local.jar MarkDuplicatesSpark -I HG00102_hGRCH38_exome_aignment.bam -O HG00102_hGRCH38_exome_aignment.marked_duplicates.bam --remove-sequencing-duplicates --spark-master local[*]

22:45:08.487 INFO  NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/mnt/c/Users/angus/Documents/Bioinformatics/GATK4/gatk-4.2.4.0/gatk-package-4.2.4.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
Dec 19, 2021 10:45:08 PM shaded.cloud_nio.com.google.auth.oauth2.ComputeEngineCredentials runningOnComputeEngine
INFO: Failed to detect whether we are running on Google Compute Engine.
22:45:08.748 INFO  MarkDuplicatesSpark - ------------------------------------------------------------
22:45:08.748 INFO  MarkDuplicatesSpark - The Genome Analysis Toolkit (GATK) v4.2.4.0
22:45:08.748 INFO  MarkDuplicatesSpark - For support and documentation go to https://software.broadinstitute.org/gatk/
22:45:08.748 INFO  MarkDuplicatesSpark - Executing as *removed* on Linux v4.19.128-microsoft-standard amd64
22:45:08.749 INFO  MarkDuplicatesSpark - Java runtime: OpenJDK 64-Bit Server VM v11.0.9.1+1-Ubuntu-0ubuntu1.18.04
22:45:08.749 INFO  MarkDuplicatesSpark - Start Date/Time: December 19, 2021 at 10:45:08 PM GMT
22:45:08.749 INFO  MarkDuplicatesSpark - ------------------------------------------------------------
22:45:08.749 INFO  MarkDuplicatesSpark - ------------------------------------------------------------
22:45:08.750 INFO  MarkDuplicatesSpark - HTSJDK Version: 2.24.1
22:45:08.750 INFO  MarkDuplicatesSpark - Picard Version: 2.25.4
22:45:08.751 INFO  MarkDuplicatesSpark - Built for Spark Version: 2.4.5
22:45:08.751 INFO  MarkDuplicatesSpark - HTSJDK Defaults.COMPRESSION_LEVEL : 2
22:45:08.751 INFO  MarkDuplicatesSpark - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
22:45:08.751 INFO  MarkDuplicatesSpark - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
22:45:08.751 INFO  MarkDuplicatesSpark - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
22:45:08.751 INFO  MarkDuplicatesSpark - Deflater: IntelDeflater
22:45:08.751 INFO  MarkDuplicatesSpark - Inflater: IntelInflater
22:45:08.751 INFO  MarkDuplicatesSpark - GCS max retries/reopens: 20
22:45:08.752 INFO  MarkDuplicatesSpark - Requester pays: disabled
22:45:08.752 INFO  MarkDuplicatesSpark - Initializing engine
22:45:08.752 INFO  MarkDuplicatesSpark - Done initializing engine
Using Spark's default log4j profile: org/apache/spark/log4j-defaults.properties
21/12/19 22:45:09 WARN Utils: Your hostname, *hostname removed* resolves to a loopback address: *IP removed*; using 172.27.147.65 instead (on interface eth0)
21/12/19 22:45:09 WARN Utils: Set SPARK_LOCAL_IP if you need to bind to another address
WARNING: An illegal reflective access operation has occurred
WARNING: Illegal reflective access by org.apache.spark.unsafe.Platform (file:/mnt/c/Users/angus/Documents/Bioinformatics/GATK4/gatk-4.2.4.0/gatk-package-4.2.4.0-local.jar) to method java.nio.Bits.unaligned()
WARNING: Please consider reporting this to the maintainers of org.apache.spark.unsafe.Platform
WARNING: Use --illegal-access=warn to enable warnings of further illegal reflective access operations
WARNING: All illegal access operations will be denied in a future release
21/12/19 22:45:10 INFO SparkContext: Running Spark version 2.4.5
21/12/19 22:45:10 WARN NativeCodeLoader: Unable to load native-hadoop library for your platform... using builtin-java classes where applicable
21/12/19 22:45:10 INFO SparkContext: Submitted application: MarkDuplicatesSpark
21/12/19 22:45:10 INFO SecurityManager: Changing view acls to: abgane
21/12/19 22:45:10 INFO SecurityManager: Changing modify acls to: abgane
21/12/19 22:45:10 INFO SecurityManager: Changing view acls groups to:
21/12/19 22:45:10 INFO SecurityManager: Changing modify acls groups to:
21/12/19 22:45:10 INFO SecurityManager: SecurityManager: authentication disabled; ui acls disabled; users  with view permissions: Set(abgane); groups with view permissions: Set(); users  with modify permissions: Set(abgane); groups with modify permissions: Set()
21/12/19 22:45:10 INFO Utils: Successfully started service 'sparkDriver' on port 35319.
21/12/19 22:45:10 INFO SparkEnv: Registering MapOutputTracker
21/12/19 22:45:11 INFO SparkEnv: Registering BlockManagerMaster
21/12/19 22:45:11 INFO BlockManagerMasterEndpoint: Using org.apache.spark.storage.DefaultTopologyMapper for getting topology information
21/12/19 22:45:11 INFO BlockManagerMasterEndpoint: BlockManagerMasterEndpoint up
21/12/19 22:45:11 INFO DiskBlockManager: Created local directory at /tmp/blockmgr-1e5e323f-dd90-4093-980c-0ceda43349c2
21/12/19 22:45:11 INFO MemoryStore: MemoryStore started with capacity 3.6 GB
21/12/19 22:45:11 INFO SparkEnv: Registering OutputCommitCoordinator
21/12/19 22:45:11 INFO Utils: Successfully started service 'SparkUI' on port 4040.
21/12/19 22:45:11 INFO SparkUI: Bound SparkUI to 0.0.0.0, and started at http://172.27.147.65:4040
21/12/19 22:45:11 INFO Executor: Starting executor ID driver on host localhost
21/12/19 22:45:11 INFO Utils: Successfully started service 'org.apache.spark.network.netty.NettyBlockTransferService' on port 33189.
21/12/19 22:45:11 INFO NettyBlockTransferService: Server created on 172.27.147.65:33189
21/12/19 22:45:11 INFO BlockManager: Using org.apache.spark.storage.RandomBlockReplicationPolicy for block replication policy
21/12/19 22:45:11 INFO BlockManagerMaster: Registering BlockManager BlockManagerId(driver, 172.27.147.65, 33189, None)
21/12/19 22:45:11 INFO BlockManagerMasterEndpoint: Registering block manager 172.27.147.65:33189 with 3.6 GB RAM, BlockManagerId(driver, 172.27.147.65, 33189, None)
21/12/19 22:45:11 INFO BlockManagerMaster: Registered BlockManager BlockManagerId(driver, 172.27.147.65, 33189, None)
21/12/19 22:45:11 INFO BlockManager: Initialized BlockManager: BlockManagerId(driver, 172.27.147.65, 33189, None)
22:45:11.804 INFO  MarkDuplicatesSpark - Spark verbosity set to INFO (see --spark-verbosity argument)
21/12/19 22:45:11 INFO GoogleHadoopFileSystemBase: GHFS version: 1.9.4-hadoop3
21/12/19 22:45:12 INFO MemoryStore: Block broadcast_0 stored as values in memory (estimated size 172.9 KB, free 3.6 GB)
21/12/19 22:45:13 INFO MemoryStore: Block broadcast_0_piece0 stored as bytes in memory (estimated size 35.5 KB, free 3.6 GB)
21/12/19 22:45:13 INFO BlockManagerInfo: Added broadcast_0_piece0 in memory on 172.27.147.65:33189 (size: 35.5 KB, free: 3.6 GB)
21/12/19 22:45:13 INFO SparkContext: Created broadcast 0 from newAPIHadoopFile at PathSplitSource.java:96
21/12/19 22:45:13 INFO MemoryStore: Block broadcast_1 stored as values in memory (estimated size 172.9 KB, free 3.6 GB)
21/12/19 22:45:13 INFO MemoryStore: Block broadcast_1_piece0 stored as bytes in memory (estimated size 35.5 KB, free 3.6 GB)
21/12/19 22:45:13 INFO BlockManagerInfo: Added broadcast_1_piece0 in memory on 172.27.147.65:33189 (size: 35.5 KB, free: 3.6 GB)
21/12/19 22:45:13 INFO SparkContext: Created broadcast 1 from newAPIHadoopFile at PathSplitSource.java:96
21/12/19 22:45:13 INFO FileInputFormat: Total input files to process : 1
21/12/19 22:45:13 INFO SparkUI: Stopped Spark web UI at http://172.27.147.65:4040
21/12/19 22:45:13 INFO MapOutputTrackerMasterEndpoint: MapOutputTrackerMasterEndpoint stopped!
21/12/19 22:45:13 INFO MemoryStore: MemoryStore cleared
21/12/19 22:45:13 INFO BlockManager: BlockManager stopped
21/12/19 22:45:13 INFO BlockManagerMaster: BlockManagerMaster stopped
21/12/19 22:45:13 INFO OutputCommitCoordinator$OutputCommitCoordinatorEndpoint: OutputCommitCoordinator stopped!
21/12/19 22:45:13 INFO SparkContext: Successfully stopped SparkContext
22:45:13.851 INFO  MarkDuplicatesSpark - Shutting down engine
[December 19, 2021 at 10:45:13 PM GMT] org.broadinstitute.hellbender.tools.spark.transforms.markduplicates.MarkDuplicatesSpark done. Elapsed time: 0.09 minutes.
Runtime.totalMemory()=419430400
java.lang.IllegalArgumentException: Unsupported class file major version 55
    at org.apache.xbean.asm6.ClassReader.<init>(ClassReader.java:166)
    at org.apache.xbean.asm6.ClassReader.<init>(ClassReader.java:148)
    at org.apache.xbean.asm6.ClassReader.<init>(ClassReader.java:136)
    at org.apache.xbean.asm6.ClassReader.<init>(ClassReader.java:237)
    at org.apache.spark.util.ClosureCleaner$.getClassReader(ClosureCleaner.scala:49)
        at org.apache.spark.util.FieldAccessFinder$$anon$3$$anonfun$visitMethodInsn$2.apply(ClosureCleaner.scala:517)
        at org.apache.spark.util.FieldAccessFinder$$anon$3$$anonfun$visitMethodInsn$2.apply(ClosureCleaner.scala:500)
        at scala.collection.TraversableLike$WithFilter$$anonfun$foreach$1.apply(TraversableLike.scala:733)
    at scala.collection.mutable.HashMap$$anon$1$$anonfun$foreach$2.apply(HashMap.scala:134)
    at scala.collection.mutable.HashMap$$anon$1$$anonfun$foreach$2.apply(HashMap.scala:134)
    at scala.collection.mutable.HashTable$class.foreachEntry(HashTable.scala:236)
        at scala.collection.mutable.HashMap.foreachEntry(HashMap.scala:40)
        at scala.collection.mutable.HashMap$$anon$1.foreach(HashMap.scala:134)
    at scala.collection.TraversableLike$WithFilter.foreach(TraversableLike.scala:732)
    at org.apache.spark.util.FieldAccessFinder$$anon$3.visitMethodInsn(ClosureCleaner.scala:500)
    at org.apache.xbean.asm6.ClassReader.readCode(ClassReader.java:2175)
    at org.apache.xbean.asm6.ClassReader.readMethod(ClassReader.java:1238)
    at org.apache.xbean.asm6.ClassReader.accept(ClassReader.java:631)
    at org.apache.xbean.asm6.ClassReader.accept(ClassReader.java:355)
    at org.apache.spark.util.ClosureCleaner$$anonfun$org$apache$spark$util$ClosureCleaner$$clean$14.apply(ClosureCleaner.scala:307)
    at org.apache.spark.util.ClosureCleaner$$anonfun$org$apache$spark$util$ClosureCleaner$$clean$14.apply(ClosureCleaner.scala:306)
        at scala.collection.immutable.List.foreach(List.scala:392)
        at org.apache.spark.util.ClosureCleaner$.org$apache$spark$util$ClosureCleaner$$clean(ClosureCleaner.scala:306)
    at org.apache.spark.util.ClosureCleaner$.clean(ClosureCleaner.scala:162)
        at org.apache.spark.SparkContext.clean(SparkContext.scala:2326)
        at org.apache.spark.SparkContext.runJob(SparkContext.scala:2100)
        at org.apache.spark.SparkContext.runJob(SparkContext.scala:2126)
        at org.apache.spark.rdd.RDD$$anonfun$collect$1.apply(RDD.scala:990)
    at org.apache.spark.rdd.RDDOperationScope$.withScope(RDDOperationScope.scala:151)
    at org.apache.spark.rdd.RDDOperationScope$.withScope(RDDOperationScope.scala:112)
    at org.apache.spark.rdd.RDD.withScope(RDD.scala:385)
    at org.apache.spark.rdd.RDD.collect(RDD.scala:989)
    at org.apache.spark.RangePartitioner$.sketch(Partitioner.scala:309)
    at org.apache.spark.RangePartitioner.<init>(Partitioner.scala:171)
    at org.apache.spark.RangePartitioner.<init>(Partitioner.scala:151)
    at org.apache.spark.rdd.OrderedRDDFunctions$$anonfun$sortByKey$1.apply(OrderedRDDFunctions.scala:62)
        at org.apache.spark.rdd.OrderedRDDFunctions$$anonfun$sortByKey$1.apply(OrderedRDDFunctions.scala:61)
        at org.apache.spark.rdd.RDDOperationScope$.withScope(RDDOperationScope.scala:151)
        at org.apache.spark.rdd.RDDOperationScope$.withScope(RDDOperationScope.scala:112)
        at org.apache.spark.rdd.RDD.withScope(RDD.scala:385)
        at org.apache.spark.rdd.OrderedRDDFunctions.sortByKey(OrderedRDDFunctions.scala:61)
        at org.apache.spark.api.java.JavaPairRDD.sortByKey(JavaPairRDD.scala:936)
        at org.broadinstitute.hellbender.utils.spark.SparkUtils.sortUsingElementsAsKeys(SparkUtils.java:165)
        at org.broadinstitute.hellbender.utils.spark.SparkUtils.sortReadsAccordingToHeader(SparkUtils.java:143)
        at org.broadinstitute.hellbender.utils.spark.SparkUtils.querynameSortReadsIfNecessary(SparkUtils.java:306)
        at org.broadinstitute.hellbender.tools.spark.transforms.markduplicates.MarkDuplicatesSpark.mark(MarkDuplicatesSpark.java:206)
        at org.broadinstitute.hellbender.tools.spark.transforms.markduplicates.MarkDuplicatesSpark.mark(MarkDuplicatesSpark.java:270)
        at org.broadinstitute.hellbender.tools.spark.transforms.markduplicates.MarkDuplicatesSpark.runTool(MarkDuplicatesSpark.java:350)
        at org.broadinstitute.hellbender.engine.spark.GATKSparkTool.runPipeline(GATKSparkTool.java:546)
        at org.broadinstitute.hellbender.engine.spark.SparkCommandLineProgram.doWork(SparkCommandLineProgram.java:31)
        at org.broadinstitute.hellbender.cmdline.CommandLineProgram.runTool(CommandLineProgram.java:140)
        at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMainPostParseArgs(CommandLineProgram.java:192)
        at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:211)
        at org.broadinstitute.hellbender.Main.runCommandLineProgram(Main.java:160)
        at org.broadinstitute.hellbender.Main.mainEntry(Main.java:203)
        at org.broadinstitute.hellbender.Main.main(Main.java:289)
21/12/19 22:45:13 INFO ShutdownHookManager: Shutdown hook called
21/12/19 22:45:13 INFO ShutdownHookManager: Deleting directory /tmp/spark-f49df209-c946-40e3-96ca-9d8125c3a733

```

As an alternative, I am trying to sort the BAM file and use MarkDuplicates instead (https://gatk.broadinstitute.org/hc/en-us/articles/4409924785691-MarkDuplicates-Picard-).This appears to work fine. I obviously want to make use of multiprocessing though so ideally would like to get the Spark implementation going.

Does anyone know what might be causing this command to fail?

Any ideas for a fix?

Any alternative software that might be useful?

Thanks for your time,

Angus

# Answer
You need to set your Java version to Java 8.

If you are using conda, run this in your conda environment:

```
conda install openjdk=8

```
> 1 votes

---
Tags: sequence-alignment, gatk, markduplicate
---


# Question
Title: How to calculate average BLOSUM62 scores?
I can understand the motive behind the BLOSUM62 matrix, this being a pairwise mutation matrix describing aggregate mutations between the 20 amino acids.

However how would you calculate the average BLOSUM62 score within a hydrophobic amino acids (K, D, E) vs hydrophilic amino acids (I, L, V)?

How about between those two groups?

# Answer
> 0 votes
I've not done this precise calculation. I would use point 3, but its complicated to implement.

Just for information BLOSUM62 is a static matrix 20x20 representing relative amino acid mutations. However, I am not exactly sure how the matrix was calculated. The algorithm looks quite simple, then a BLOSUM62 score can be obtained by comparing the observed matrix to the static matrix, this would provide a score. The paper is here.

Maximum likelihood approaches require a phylogeny.

There are three possibilities I can think of, the first one wouldn't work under a maximum likelihood approach:

1. Remove all hydrophobic residues from the alignment, check the BLOSUM62 score, remove all hydrophilic check the score etc... again a maximum likelihood approach is not meaningful in this example, but a BLOSUM62 score might be okay;
2. Recode the amino acids to generate 6-state amino acid data, or recode hydrophobics and in a separate data set hydrophilics and leave the rest uncoded. I've never done it but at a guess `I, L, V` could be represented by just one amino acid, e.g. `L`. If the observed BLOSUM62 matrix algorithm was calculated on the observed data the new amino acid e.g. `L` would have an `average` score.
3. Alternatively, manipulate the BLOSUM62 matrix, however its no longer a BLOSUM62 matrix, but it is highly amenable to maximum likelihood analysis. Thus I \<-\> V and I \<-\> L is set to zero (or massively negatively weight them) check likelihood, then hydrophilics are set to zero check likelihood. Do the same with the observed versus BLOSUM62 matrix as the control.

Point 3 what a maximum likelihood approach would do is say how close a given matrix is to observed data set for a given tree under a likelihood probability, the smallest is the best. So if you've got loads of matrices its a really cool approach. There is a membrane protein matrix - which is likely heavily skewed towards hydrophilicity. Basically, I'd test all the matrices - but not recode or delete the data set in any way. The matrix that is closer to the data set using likelihood is the best representation of its hydrophilicity vs. hydrophobicity.

**Note** Just finally mention you can simply calculate the observed amino acid matrix for that data set (via a program) and then average the amino acids. That might be easier and using parsimony its very quick but you'd need a tree. The underlying metric (or metric behind the metric) in all cases is mutations.

---
Tags: proteins, sequence-alignment, mutations
---


# Question
Title: Plink 1.9 equivalent of --score sum in Plink 2.0
I have previously performing Allelic scoring in Plink 1.9 using:

```
plink \
--score sum sumstats.txt\

```

However, it seems this has changed in Plink 2.0. What is the equivalent function for Allelic scoring?

# Answer
> 3 votes
```
plink2 --score cols=+scoresums

```

See https://www.cog-genomics.org/plink/2.0/general\_usage#colset and https://www.cog-genomics.org/plink/2.0/score for details.

# Answer
> 1 votes
I don't have enough reputation for a comment but the place where the option to add scoresums is listed is here https://www.cog-genomics.org/plink/2.0/formats#sscore

---
Tags: plink
---


# Question
Title: Gencode PolyA feature annotation GTF gene_id not the same as Comprehensive gene annotation GTF
I have downloaded the Gencode PolyA feature annotation GTF here. The first 10 lines are like this:

```
##description: evidence-based annotation of the human genome (GRCh38), version 41 (Ensembl 107) - polyA features
##provider: GENCODE
##contact: gencode-help@ebi.ac.uk
##format: gtf
##date: 2022-05-12
chr1    HAVANA  pseudo_polyA    134957  134962  .   -   .   gene_id "383460"; transcript_id "383460"; gene_type "pseudo_polya"; gene_name "383460"; transcript_type "pseudo_polya"; transcript_name "383460"; level 2;
chr1    HAVANA  polyA_site  141474  141475  .   -   .   gene_id "383461"; transcript_id "383461"; gene_type "polya_site"; gene_name "383461"; transcript_type "polya_site"; transcript_name "383461"; level 2;
chr1    HAVANA  polyA_signal    141477  141482  .   -   .   gene_id "383462"; transcript_id "383462"; gene_type "polya_signal"; gene_name "383462"; transcript_type "polya_signal"; transcript_name "383462"; level 2;
chr1    HAVANA  pseudo_polyA    184934  184939  .   -   .   gene_id "383463"; transcript_id "383463"; gene_type "pseudo_polya"; gene_name "383463"; transcript_type "pseudo_polya"; transcript_name "383463"; level 2;
chr1    HAVANA  pseudo_polyA    258530  258535  .   -   .   gene_id "383464"; transcript_id "383464"; gene_type "pseudo_polya"; gene_name "383464"; transcript_type "pseudo_polya"; transcript_name "383464"; level 2;

```

In this GTF file, what does the `gene_id` mean? In the whole genome GTF (same build) the `gene_id` is clearly the ENSEMBL gene\_id:

```
##description: evidence-based annotation of the human genome (GRCh38), version 41 (Ensembl 107)
##provider: GENCODE
##contact: gencode-help@ebi.ac.uk
##format: gtf
##date: 2022-05-12
chr1    HAVANA  gene    11869   14409   .   +   .   gene_id "ENSG00000223972.5"; gene_type "transcribed_unprocessed_pseudogene"; gene_name "DDX11L1"; level 2; hgnc_id "HGNC:37102"; havana_gene "OTTHUMG00000000961.2";
chr1    HAVANA  transcript  11869   14409   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_name "DDX11L1"; transcript_type "processed_transcript"; transcript_name "DDX11L1-202"; level 2; transcript_support_level "1"; hgnc_id "HGNC:37102"; tag "basic"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000362751.1";
chr1    HAVANA  exon    11869   12227   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_name "DDX11L1"; transcript_type "processed_transcript"; transcript_name "DDX11L1-202"; exon_number 1; exon_id "ENSE00002234944.1"; level 2; transcript_support_level "1"; hgnc_id "HGNC:37102"; tag "basic"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000362751.1";
chr1    HAVANA  exon    12613   12721   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_name "DDX11L1"; transcript_type "processed_transcript"; transcript_name "DDX11L1-202"; exon_number 2; exon_id "ENSE00003582793.1"; level 2; transcript_support_level "1"; hgnc_id "HGNC:37102"; tag "basic"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000362751.1";
chr1    HAVANA  exon    13221   14409   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_name "DDX11L1"; transcript_type "processed_transcript"; transcript_name "DDX11L1-202"; exon_number 3; exon_id "ENSE00002312635.1"; level 2; transcript_support_level "1"; hgnc_id "HGNC:37102"; tag "basic"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000362751.1";

```

What I want to do eventually is to assign the PolyA sites to ENSEMBL gene\_ids.

# Answer
According to the GENCODE ftp README file, the value of the 'gene\_id', 'transcript\_id', 'gene\_name' and 'transcript name' fields corresponds to a random identifier. The polyA features are not directly associated to any gene or transcript when manually annotated by Havana.

There is other file in the GENCODE release that includes links between polyA features and transcripts: https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode\_human/release\_41/gencode.v41.metadata.PolyA\_feature.gz

```
ENST00000485748.5       3232    3233    chr1    1496201 1496202 +       polya_site
ENST00000485748.5       3211    3216    chr1    1496180 1496185 +       polya_signal
ENST00000474481.1       3367    3368    chr1    1496201 1496202 +       polya_site
ENST00000474481.1       3346    3351    chr1    1496180 1496185 +       polya_signal
ENST00000308647.8       2114    2115    chr1    1496201 1496202 +       polya_site
ENST00000308647.8       2093    2098    chr1    1496180 1496185 +       polya_signal

```

```
22. gencode.vX.metadata.PolyA_feature.gz:
  Manually annotated polyA feature overlapping the transcript 3'-end.
  1 - transcript id
  2 - transcript-based start coordinate of the polyA feature
  3 - transcript-based end coordinate of the polyA feature
  4 - polyA feature chromosome
  5 - polyA feature start coordinate
  6 - polyA feature end coordinate
  7 - polyA feature strand
  8 - polyA feature type ("polyA_site", "polyA_signal", "pseudo_polyA")

```
> 2 votes

---
Tags: sequence-annotation, ensembl, gtf, gencode
---


# Question
Title: Coordinate numbering between forward (+) and reverse (-) strand in Tombo
My question is related to the output of Tombo re-squiggling algorithm on nanopore data. Given the position/coordinate on one strand (+ or -), I want to find the coordinate of its exact location on the opposite strand. For example, let's say I want to find data/signal on specific coordinates on both strands on Tombo path after re-squiggling. Let's say position 13936 - 14036 on genome on + and - negative strands. I wrote the following code but the sequence obtained from Tombo is not similar at all. I want to know how I can find the coordinate corresponding to a given position on the opposite strand.

```
# one read from the + strand
mapped_start, mapped_end, read_id, fast5_path, _, _ = reads_per_win_pos.loc[rd_id, :]
>>> start
13936
>>> end
13942
>>> mapped_start, mapped_end, read_id, fast5_path, _, _ = reads_per_win_pos.loc[rd_id, :]
>>> mapped_start
11205
>>> mapped_end
20084
>>> f = h5py.File(fast5_path, 'r')
>>> tombo_events = list(f['Analyses/RawGenomeCorrected_000/BaseCalled_template/Events'])
>>> forward_seq = ''.join([e[4].decode() for e in tombo_events])
>>> start_idx = start - mapped_start
>>> start_idx
2731
>>> win_seq = seq[start_idx: start_idx + 100]
>>> forward_win_seq = forward_seq[start_idx: start_idx + 100]
>>> # one other overlapping read on the other strand (-)
...
>>> mapped_start, mapped_end, read_id, fast5_path, _, _ = reads_per_win_neg.loc[rd_id, :]
>>> mapped_start
10979
>>> mapped_end
20091
>>> f = h5py.File(fast5_path, 'r')
>>> tombo_events = list(f['Analyses/RawGenomeCorrected_000/BaseCalled_template/Events'])
>>> reverse_seq = ''.join([e[4].decode() for e in tombo_events])
>>> start_idx = start - mapped_start
>>> start_idx
2957
>>> reverse_win_seq = reverse_seq[start_idx: start_idx + 100]

```

Now the forward\_win\_seq and reverse (or even the reverse complement) of it are not similar:

```
>>> forward_win_seq
'GAGACGGAGCAGACCCATCTGCTACTGCCCTTTCTATAATAACTAAAGTTAGCTGCCCTGGACTATTCACCCCCTAGTCTCAATTTAAAAAGATCCCCAT'
>>> reverse_win_seq
'GGGCCAGGTGCCTGAGATCCATGTTCCATCCTACCTGCCTGACCTGCCCGGCATTGCCAACGACCTCATGTACAGTGCCGACCTGGGCCCCGGCATTGCC'
>>> complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}
>>> com_seq = ''.join([complement[base] if base in complement.keys() else base for base in reverse_win_seq])
>>> com_seq[::-1]
'GGCAATGCCGGGGCCCAGGTCGGCACTGTACATGAGGTCGTTGGCAATGCCGGGCAGGTCAGGCAGGTAGGATGGAACATGGATCTCAGGCACCTGGCCC'

```

# Answer
> 2 votes
When **Tombo re-squiggle** and aligns a read it adds two attributes to it: **mapped\_start** and **mapped\_end**. Tombo also gives chromosomes and strands.

Imagine you are looking for **some value of a signal\[s\] in position x to y on a chromosome on both strands.** In this range, you probably have reads that span the positive strand (forward) and the reads that cover the negative strand (reverse).

Now let a given read that has mapped to a positive strand and covers the entire \[x, y\], have attributes: **pos\_mapped\_start** and **pos\_mapped\_end**. Similarly let the reads mapped to negative strands covering \[x, y\] have attributes: **neg\_mapped\_start** and **neg\_mapped\_end**.

For reads mapped to + strand:

```
start_idx = x- pos_mapped_start
end_idx = y - pos_mapped_start + 1
required_signal= read_signal[start_idx: end_idx]

```

For reads mapped to - strand:

```
start_idx = neg_mapped_end - y- 1
end_idx = neg_mapped_end - x
required_signal= read_signal[start_idx: end_idx]

```
* Note that x and y are **inclusive** and the length of the required signal = y - x + 1.

---
Tags: nanopore, read-mapping
---


title: need bam file for pilon

text:
# Question
Title: need bam file for pilon
I just ran an assembly on yeast genomes using Flye and I want to polish those assemblies with Pilon but it requires a sorted BAM file.

How do I make a BAM file of the resulting assembled.fasta?

# Answer
Run a short-read mapper. For example:

```
bwa index assembled.fasta
bwa mem -pt16 assembled.fasta read1.fq.gz read2.fq.gz \
  | samtools sort -m4G -@4 -o align.bam -

```
> 2 votes

---
Tags: sam, sequence-analysis, assembly, conversion, genome-index
---


# Question
Title: How to chop fasta / bed /peak file genomic segments into smaller fixed or custom genomic intervals
I am trying to split my *.bed/*.fa file genomic segments into arbitrary smaller overlapping intervals:

Consider a typical line/row in my \*.bed file as follows:

chrom. $\ \ \ \ \ $ start $\ \ \ \ \ \ \ \ \ \ \ \ $ end

chr7 $\ \ \ $ 127471196 $\ \ \ $ 127472363$\ \ \ $ ...

The start and end genomic positions include 1167 base pairs for this line. I need to extend this line, for example, into multiple lines where each line covers exactly 500bp like below:

chr7 $\ \ \ $ 127471196 $\ \ \ $ 127471696 $\ \ \ $ ...

chr7 $\ \ \ $ 127471863 $\ \ \ $ 127472363 $\ \ \ $ ...

While each line covers 500pb, between the two new lines they both include the same "range" of genomic positions as the original line including some overlap which is fine (and desired)

Is there an existing tool (perhaps part of bedtools?) that can do something like this?

I know there is the **bedops** tool that does ALMOST this using the --chop option, where it takes the original segment x1:x2 and chops the original segment into specified fixed bp length intervals, however, if there is any left over towards the end is creates a new non-standard length segment (i.e. x1:x2=10bp and chop length=3, it will give you 4 segments that are length 3, 3, 3 and 1 - what I am looking for is something that adjusts for the last remainder and outputs segments of length 3,3,3 and 3 where the last length 3 segment is right adjusted, while the first three segments are left adjusted or better yet it gives the user the choice of precise redefinition of segment coordinates y1,y2,y3 and y4). Thank you.

# Answer
You can run `--chop` with `-x` to excise the last piece, adding the `--stagger` to size the chopped pieces as needed, e.g.:

```
$ echo -e 'chr1\t0\t10' | bedops --chop 3 --stagger 2 -x -
chr1    0       3
chr1    2       5
chr1    4       7
chr1    6       9

```

Adjust parameters to taste.

https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html#chop-w-chop

https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html#stagger-stagger

https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html#exclude-x
> 2 votes

---
Tags: python, fasta, sequence-alignment, chip-seq, bedtools
---


# Question
Title: Downloading genomic protein files from accessions in Python
I am trying to download the `_protein.faa.gz` files for genomes given their accession numbers through Python. Ideally, I would like to do this without third party libraries. Essentially what I have is a list of just the GCA or GCF accessions. The issue with the ftp site is that it includes the project name along with the accession, which I do not have ahead of time. It would be perfect if I could run downloads from within Python that matched a pattern something like this:

`https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/145/295/GCF_000145295.1_*/*_protein.faa.gz`

Is it at all possible to do this? Or is it possible to run something similar *via* efetch in Biopython?

Thanks in advance!

# Answer
Concept code only

```
from Bio import Entrez

Entrez.email = "m@M__"
handle = Entrez.efetch(db="protein", id="GCF_000145295", "GCF_000145294", "GCF_000145293", retmode="text",rettype="gb") # check the quotes
records = Entrez.parse(handle)
handle.close()
for record in records:
    print (record)

```

You can write a loop, flatten it and then insert it to the `id=` tag.

```
mylist = ['GCF_000145295', 'GCF_000145294', 'GCF_000145293']
flattened = ','.join(mylist)

```
> 1 votes

# Answer
The above answer is right - here are some more links that describe what you want to do in more detail:

You're going to want to use the BioPython package as it has functions that do exactly what you're looking for.
> 3 votes

---
Tags: python, phylogenetics
---


# Question
Title: Building a group representation of genes' ontologies
I have a list of genes (about 2000) from lines of **drosophila subobscura**, in the following format:

* LOC117900589
* LOC11788959
* LOC11790331

I'm trying to represent this set of genes in a useful way.

So, I was recommended to represent the ontology function of these genes in a graphical way. I was recommended this tool: http://cbl-gorilla.cs.technion.ac.il/

One inputs a list of genes and it outputs this kind of graph:

The problem is that, this tool, GOrilla only accepts genes from **drosophila melanogaster**.

I have to either:

* **a)** convert the genes id to a homologous gene in **melanogaster**
* **b)** find an alternative tool

A solution to either a) or b) would be helpful. Or an alternative solution.

For a) I tried:

* Ncbi's **homologene**: which doesn't seem to be working for a the genes ID I tried so far (I didn't try them all)
* Blasting the gene sequence on melanogaster's genome and looking through the list manually to find an homologous gene.

So far these methods have been impractical for the number of genes I have.

# Answer
> 1 votes
Have you tried DAVID for functional annotation (Gene Ontology and KEGG terms)?

I usually use DAVID to extract GO and KEGG terms and plot the results using R.

If you also want to plot a protein-protein interaction (PPI) network for your gene list you can use STRING database to extract interactions and network and use Cytoscape for visualizing this result.

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


# Question
Title: Merge/Join one data frame (X) to multiple data frames separately
Let us say that we have some data frames named S1,S2,S3... that I want to join independently by column names to dataframe X. This is how they look like:

```
S1 <- data.frame(
  values = c (1:5),
  numbers = c(11:15))
S1

S2 <- data.frame(
  values = c (1:5),
  numbers = c(11:15))
S2

S3 <- data.frame(
  values = c (1:5),
  numbers = c(16:20))
S3

```

#Dataframe X

```
X <- data.frame(
  values = c (26:30),
  numbers = c(501:505))
X

S1
S2
S3
X

```

Using rbind, it is really easy to merge/join them according to column names. As below:

```
S1_X<- rbind(S1,X)

S1_X

S2_X<- rbind(S2,X)

S2_X

S3_X<- rbind(S3,X)

S3_X

```

Now the problem is that I have 500 plus data frames (S1,S2,S3,S4,S5.....S500) that I want to join to dataframe "X" by any function. For example rbind. Is it possible to use some loop/function for that.

What I have best tried is create a list (nested list of dataframes) that I can later split into individual dataframes:

```
output = lapply(list(S1=S1,S2=S2,S3=S3),rbind,X)

```

The output is what I expected. So it works The only problem is, how do I specify S1 to S500?

Thank you in advance!

# Answer
> 1 votes
You can use For Loop here. Here I have created code using your algorithm, first created a list of merged dataframes and then split the list into separate dataframes.

```
S1 <- data.frame(values = c (1:5), numbers = c(11:15))
S2 <- data.frame(values = c (1:5), numbers = c(11:15))
S3 <- data.frame(values = c (1:5), numbers = c(16:20))
#Dataframe X
X <- data.frame(values = c (26:30), numbers = c(501:505))

#defining a NULL object
S=NULL
#creating a list of merged dataframes
for (i in 1:3){
S[i] = list(rbind(get(paste("S",i,sep='')),X))
names(S)[i] = paste("S",i,"_X",sep='')}

#extracting dataframes
for(i in 1:3){
assign(paste("S",i,"_X",sep = ""),S[[i]])
}

```

---
Tags: r, dataframe
---


# Question
Title: WCGNA - Relate modules with Y features when the % of variance explained of each eigengen is low
I'm doing a WCGNA analysis (signed network) on microbiome 16S data. I have transformed counts to centeres log-ratio transformed data (CLR) to address the compositional characteristics of the data and have obtained a pretty decent cluster (around 7 modules, 20 taxa each).

I correlated the different modules eigengenes to my features of interest using both bicor and pearson correlation measures. Overall, the results I'm getting make sense but I'm having some questions regarding the method: The % of variance explained of each eigengen is around 30%. Is this value normal? Is it ok that I'm using only the PC1 of the module as a summary of the module itself when variances explained of the PC1 are so low?

Thx!!

# Answer
**Part 1** Just to clarify on the specific issue of using PC1 alone given the reported result ...

Without knowing the underlying ecology, which is important, 30% is too low to any draw conclusions in PCA-based analysis, particularly because this could be a really complex niche, ecotope, ecoregion, etc ...

Thus this specific result from the unsupervised learning shouldn't be used as the basis for any downstream analysis. This is because the remaining 70% of the variance could overturn PC1, where PC1 could be separating a minor feature of the community ecology.

Thus, its informative because it states the first component is insufficient to proceed - however I stress its difficult to give a clear answer without knowing the expectedly complexity of the meta-genomic sample.

I suspect its a nice result, but in broad general terms its not a good idea.

---

**Part 2**

> of each eigengen is around 30%.

If PC1-PC3 is 90%, then thats good and maybe its not so complex. Personally, I'd be delighted with that result. However, the inference is that the clustering is not clear. I would simply suggest further data mining is required.
> 1 votes

# Answer
Proportion of variance explained of about 30% is indeed on the low side, but nothing out of the ordinary. I have no experience with 16S data; in gene expression data, between 30 and 50% is actually considered good. PVE byt the eigengene depends strongly on the number of samples; the more samples, the lower the PVE will be.

You can try to calculate additional singular vectors and their PVE. Usually this is much lower than the eigengene (1st singular vector). So I would focus more on the biological significance - do the modules make biological sense?
> 0 votes

---
Tags: rna-seq, pca, wgcna, microbiome, 16rrna
---


# Question
Title: Visualize RNA co-folded secondary structure with Python
I would like to visualize the secondary structure of co-folded RNA strands in Python, preferably in a Jupyter notebook. What are the recommended tools? Can anyone provide an example?

# Answer
> 1 votes
... . I know both aspects of this question and I'm not aware of an interactive Python package using e.g. iPython magic commands that will do this. The thing you could try is `draw_RNA` https://github.com/DasLab/draw\_rna

`draw_RNA` will operate on Jupyter, but I've not tried it.

A very, very good package is VARNA (its Java BTW). The `VARNA` group make further recommendations from there for complex structures. I would vouch for this package for complex structures.

A good approach is to output in Stockholm 1.0 format. It then just drops and drags into `VARNA`. `Dot-bracket` format also works and I would assume `ct` format works too.

A good reason why in-line graph imaging isn't the best solution for complex structures, i.e. the in-line graphs in a Python debugger or Jupiter is complex RNA structures need to be interactive because you need to move them: they "clump" together and have to manually separate them.

As a final note Forna may be of benefit,

---
Tags: python, visualization, rna, secondary-structure, jupyter
---


# Question
Title: How can I work out gram staining information from the species name?
I have a very long list of prokaryote organism names. I would like to know their **gram staining** (gram positive or gram negative) information.

Please let me know if there is any database that keeps this information.

I am also trying to find a program that can find gram staining information via organism name only as my organism list is long.

Edit: for each species the NCBI taxonomy lineage is also available

# Answer
The GOLD database can be used to retrieve a list of species names according to their Gram staining: go to the "Advanced search" page and select the either "gram-" or "gram+" under "Organism Fields" -\> "Gram Stain". Then on the results page click on the number under "Organisms". On the results page you can select to have 1000 results per page, meaning that you would need to do some manual work.

I should add that as far as I know this information is based on what people submit when they register their genome projects, so I wouldn't expect this list to be curated or 100% accurate.

If you have access to the genome of the species you are interested in (you could download a representative strain for each species using ncbi-genome-download), you could use traitar, which can predict the gram staining, among other interesting broad phenotypes.
> 3 votes

# Answer
AMR R package can assign gram staining based on taxonomy. See `mo_gramstain(x, language = get_AMR_locale(), ...)`.
> 1 votes

---
Tags: database, public-databases, microbiology
---


# Question
Title: KEGG retrieve hierarchy for organism in biopython
I want to obtain a direct graph from KEGG as well as the genes present in a single term. While it is pretty easy to get list of genes/compounds that are part of a specific category, I am kinda stucked to obtain the hierarchical relationships between terms. The docs on biopython.KGML are not really helpful. This is what I have so far

```
from Bio import SeqIO
import Bio.KEGG.REST as kegg
from Bio.KEGG.KGML import KGML_parser
from Bio.Graphics.KGML_vis import KGMLCanvas
from Bio.Graphics.ColorSpiral import ColorSpiral
from itertools import chain

k =  kegg.kegg_list('pathway', 'pae').read()
todf = []
for x in k.split('\n'):
    if x=='':
        next
    dd = x.split(' - ')
    dd = [v.split('\t') if '\t' in v else v for v in dd]
    dd = list(chain.from_iterable(dd))[:2]
    try:
        ll= kegg.kegg_get(dd[0]).read().split('GENE')[-1]
        ll = [x.lstrip() for x in ll.split('COMPOUND')[0].split('\n')]
        ll= [x.rstrip().split('\t') for x in ll]
        for gn in ll:
            gn.extend(dd)
        todf.append(pd.DataFrame(ll))
    except IndexError:
        pass
todf = pd.concat(todf)
# remove compounds and random strings
todf= todf[todf[0].str.startswith('PA')]
todf[0] = [x.split(';')[0] if ';' in x else x for x in todf[0]]
todf[0] = todf[0].str.split('  ')
todf=todf.explode(0)

```

and now I have the protein to term relationship but not sure how to get the hierarchy between terms. Any idea?

# Answer
Ok realized it was much easier than expected. On KEGG it is possible to DL a json file, which is intrinsically ordered. So this parses it into a df

```
import pandas as pd
import json

todf = []

kegg = json.load(open('meta/KEGG_pae00001.json'))
for main in kegg['children']:
    for broad in main['children']:
        for sub in broad['children']:
            for gn in sub.get('children', [None]):
                if gn:
                     nm = gn['name'].split(';')[0].split(' ')
                     todf.append([nm[0], nm[1], sub['name'], broad['name'], main['name']])
                else:
                     # some terms have no children
                     pass
                     #print(sub)

todf = pd.DataFrame(todf)
todf.columns=['Locus', 'GeneName', 'Cat1', 'Cat2', 'Cat3']
todf.to_csv('meta/kegg_parsed.csv')

```
> 2 votes

---
Tags: python, kegg
---


# Question
Title: Parsing a genbank file and outputting specific feature information to a csv using BioPython
So I am trying to parse through a genbank file, extract particular feature information and output that information to a csv file. The example genbank file looks like this:

SBxxxxxx.LargeContigs.gbk

```

LOCUS       scaffold_31            38809 bp    DNA              UNK 01-JAN-1980
DEFINITION  scaffold_31.
ACCESSION   scaffold_31
VERSION     scaffold_31
KEYWORDS    .
SOURCE      .
      ORGANISM  .
COMMENT     ##antiSMASH-Data-START##
            Version      :: 6.1.1
            Run date     :: 2022-09-21 11:09:55
            ##antiSMASH-Data-END##
FEATURES             Location/Qualifiers
            protocluster    26198..38809
                            /aStool="rule-based-clusters"
                            /category="terpene"
                            /contig_edge="True"
                            /core_location="[36197:37079](-)"
                            /cutoff="20000"
                            /detection_rule="(Terpene_synth or Terpene_synth_C or
                            phytoene_synt or Lycopene_cycl or terpene_cyclase or NapT7
                            or fung_ggpps or fung_ggpps2 or trichodiene_synth or TRI5)"
                            /neighbourhood="10000"
                            /product="terpene"
                            /protocluster_number="1"
                            /tool="antismash"

```

Now for the output file, I want to create a csv with 3 columns. One column will have the Scaffold information (ie. scaffold\_31), the second column will have the category value in the protocluster feature (ie. /category = "terpene") and the third column will have the product value in the protocluster feature (ie. /product="terpene")

This is what I have so far for code. I know I can sort through the feature.qualifiers in the protocluster feature to get the category and product. I am not sure how to extract the scaffold information. I am completely new to parsing through gene bank files so have little knowledge in this domain. Thanks in advance for any assitance!

import Bio
from Bio import SeqIO
import os

input = "/Path to SBxxxxxx.LargeContigs.gbk"
output = open("output.csv", "w")

if not os.path.exists(output):
     for record in SeqIO.parse(input, "genbank")
          for feature in record.features:
              if feature.type == "protocluster" and "category" and "product" in feature.qualifiers:
                  outfile = feature.qualifiers["category"][0] + "," + feature.qualifiers["product"][0] + "\n"
                  output.write(outfile)

```
</code></pre>

# Answer
> 2 votes
Depending on which field you want to pull the "scaffold\_31" text from, you have a few options:

| GenBank Text | BioPython object |
| --- | --- |
| `LOCUS` | `record.name` |
| `DEFINITION` | `record.description` |
| `ACCESSION` | `record.accessions[0]` |
| `VERSION` | `record.version` |

Python's built in `dir()` function is handy for figuring out this kind of thing.

If you're working with a draft flat file (like BankIt gives you just before submitting) note that some of those are placeholders that get updated with the actual accession info when it's finalized. So your "scaffold\_31" text will only show up I think in the DEFINITION line in the end if I remember right.

# Answer
> 0 votes
I've used SARS-CoV-2 (Genbank: PA544053), because there was no Genbank entry given in the OPs question. If my example is representative (might not be) I think its about the object attributes. To understand the object I listed its attributes,

```
print(record[0].__dict__.keys())

```

Output

> dict\_keys(\['\_seq', 'id', 'name', 'description', 'dbxrefs', 'annotations', '\_per\_letter\_annotations', 'features'\])

The docs and @jesse's very kind response says there's a 'accession' attribute (Biopython docs below). @Jesse did mention `dir()` which was cool. I couldn't find `record[0].accession` or perhaps `record[0].accessions` and the OP might have had the same problem.

In my example *there is an `'annotations'` attribute* and beneath that was `'accession'` accessed via,

```
print(record[0].annotations['accessions'])

```

Output

> \['PA544053'\]

There it is ... Here's the full code including the CSV package, I'm using `efetch` so it'll just copy and paste and run.

```
from Bio import Entrez
from Bio import SeqIO
import csv
from pathlib import Path

fileout = Path('/mypath/mydir', 'test.csv')
handle = Entrez.efetch(db="nucleotide", id="PA544053", rettype="gb",retmode="text") #,
record = [rec for rec in SeqIO.parse(handle,'genbank')]
fieldnames = ['Accession', 'db_xref', 'organism']

# print(record[0].__dict__.keys())
print(record[0].annotations['accessions'])

with open (fileout, 'w') as f:
    csvwriter = csv.writer(f, delimiter="\t", quotechar='"')
    csvwriter.writerow(fieldnames)
    for feature in record[0].features:
        if feature.type == 'organism' or "mol_type" and "db_xref" in feature.qualifiers:
            csvwriter.writerow([''.join(record[0].annotations['accessions']), ''.join(feature.qualifiers['db_xref']), ''.join(feature.qualifiers['organism'])])

```

Output:

```
Accession      db_xref       organism
PA544053       taxon:2697049 Severe acute respiratory syndrome coronavirus 2

```

---

**Biopython docs** A likely reason for the question is the missing attribute is described in the official docs. Using a `GenBank` object (not `SeqIO`) there is certainly an `accession` attribute, https://biopython.org/docs/1.75/api/Bio.GenBank.html

The official Biopython page ...

> ```
> from Bio import GenBank
> with open("GenBank/NC_000932.gb") as handle:
>     record = GenBank.read(handle)
>     print(record.accession)
>
> ```

I think the basis of the question is to associate the accession number with the biochemical/genetic info. Note, I don't know the difference between `SeqIO` and `GenBank` objects.

---
Tags: python, genbank, parsing
---


title: convert GT:PL format to GT:GP format

text:
# Question
Title: convert GT:PL format to GT:GP format
I have a vcf info formatted in GT:PL.

`chr22 49994037 . G A 345.64 PASS AC=1;AF=0.500;AN=2;BaseQRankSum=-3.605;DP=35;ExcessHet=3.0103;FS=1.377;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=10.17;ReadPosRankSum=-0.035;SOR=0.392 GT:AD:DP:GQ:PL 0/1:20,14:34:99:353,0,588`

I want to convert it into GT:GP. But I have to write a script for that. Is there any way to do it fast and simple?

Note:

GT:Genotype

PL: Phred-scaled Genotype Likelihoods

GP: Genotype Probabilities

check the link below for all the info tags

https://samtools.github.io/hts-specs/VCFv4.2.pdf

# Answer
Currently you have **phred-scaled genotype likelihoods** or **PL** values. Genotype likelihoods give the likelihood of observing the read data (i.e. the number of reads which map to either the ref or alt allele), given a particular underlying genotype (0/0, 0/1, 1/1 at a bi-allelic locus); $GL = p(G|D)$. Phred scaling simply scales it like $PL = -10 * log(GL)$. These are usually given by a genotype caller.

**Genotype probabilities** are slightly different in that they are posterior probabilities that the genotype is correct, given either other genotypes from a population or trio and a prior distribution for the allele frequencies. Therefore, you can't trivially go from likelihoods to probabilities without accounting for this external information and the prior. As suggested by JRodrigoF, you should use something like CalculateGenotypePosteriors to do this, maybe including 1000 genomes allele counts.
> 1 votes

---
Tags: vcf, format-conversion, conversion
---


title: Searching for HLA-B in DNA results

text:
# Question
Title: Searching for HLA-B in DNA results
I'm trying to find the HLA-B\*15:01 variant in my DNA results, prompted from this research paper: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8142661/#:~:text=HLA%2DB\*15%3A01%20is%20strongly%20associated%20with%20asymptomatic,mechanism%20underlying%20early%20viral%20clearance.

I have DNA test results from Dante Labs that contain the following files

```
"*.bam"
"*.bam.bai"
"*.cnv.vcf.gz"
"*.cnv.vcf.gz.tbi"
"*.filtered.indel.vcf.gz"
"*.filtered.indel.vcf.gz.tbi"
"*.filtered.snp.vcf.gz"
"*.filtered.snp.vcf.gz.tbi"
"*1_SA_L001_R1_001.fastq.gz"
"*SA_L001_R2_001.fastq.gz"

```

I tried loading this on the Genome Explorer on Sequencing.com, but it couldn't find HLA-B at all, only HLA-A and HLA-C.

Then I tried loading the VCF files in IGV. I was able to filter to Chromosome 6 and HLA-B, which resolved to position

```
chr6:31,319,652-31,326,956

```

The CNV VCF says no variables found for chr6:31,319,652-31,326,956 The filtered Indel VCF did produce some variants for this position, but I have no clue how to search deeper specifically for HLA-B\*15:01. Here is a screenshot:

The filtered SNP VCF produced *alot* more variants for this position, but again, I don't know how to go deeper than this.

I'm not really sure what to do from here?

I don't have data on what the reference was, so that is unknown to me. I will reach out to Dante Labs about that, is there some sort of standard that is supposed to be used, or is it different for every genome provider?

How do I actually perform the check in IGV assuming the reference did not have it and the VCF did. I also assume this VCF does not have this allele, but I frankly am lost into how to confirm it to the "last mile" through the VCF so to speak, or even just go straight to BAM

# Answer
Typical SNP calling wouldn't give you good results on classical HLA genes. Ignore the VCF. Try a dedicated HLA genotyper on raw reads instead. I heard optitype is a good one for exome or whole-genome data. I haven't used it myself, though.
> 2 votes

# Answer
The specific question and the easiest question is the absence of HLA-B\*15:01 from your sample and that is not unexpected at all.

The maximum frequency of this allele is 40% found in "China Beijing pop 2" , this however is really unusually high and was published by Guang Yang et al (2009, IJI 37, 205-212).

Mostly this allele is present at 5% for moderate frequency of occurrence for example Europe, such as in Germany here. Moreover, in some populations the frequency is very low \<\<1%, e.g. India or Israel.

So a sample genome at random the will not have the B\*15:01.

---

The more complicated issue is the VCF and the reference genome. This is because the basis on whether an allele is present or absent is a comparison against the reference genome which will generate the VCF. Thus understanding the HLA-B alleles of the reference genome is really important. You need a reference genome which lacks HLA-B\*15:01, therefore if your VCF also lacks the allele, its missing. However, if HLA-B\*15:01 is present in the reference genome your VCF will not present the allele if its present.

The probability is that the allele is missing from both the reference and sample genome.

---

**Comments**

**Easy way** The easiest way is search the IGV for LT618832.1 (HLA-B\*15:01 reference), IGV has a motif finder, here https://software.broadinstitute.org/software/igv/motif\_finder

**Hard way** This is how I'd do it ..

* convert the bam to fasta then
* run a blast against Genbank reference LT618832.1

I'd section the bam file into the locus of interest (e.g. chromosome 6) e.g. ... `samtools view ... ChrName: ... > chromo6.bam`: you'd need to cross check how to do this last bit ... then via `samtools >1.3`,

```
samtools fasta chromosome6.bam > chromosome6.fasta

```

Finally you download a standalone version of `blast`, I assume its on `-c bioconda` ...

```
blastn -query HLAB1501.faa -db chromosome6.fa -out myhits.txt

```

From memory I don't think you need to `makeblastdb` for a fasta onto a fasta db.

You'd need to map the SNPs in introns parsing (the coordinates are in the Genbank reference), or optically against the Genbank reference onto the `blast` output.

You may want to cross check the HLA-B\*15:01 against https://www.ncbi.nlm.nih.gov/nuccore/?term=HLA-B%2A15%3A01

\[but be careful about the introns\]

`motif_finder` seems easier and accommodates potential heterozygosity.
> 1 votes

---
Tags: gene, genome, igv
---


# Question
Title: select highly variable genes out of dataframe
lets say I have the following dataframe:

```
Ensembl_ID          A   B   C   D   E   F   G
ENSG00000236601.2   0   0   0   1   0   0   1
ENSG00000237094.12  2   187 442 109 88  144 486
ENSG00000269732.1   0   0   0   0   0   0   0
ENSG00000284733.1   0   0   0   0   0   0   0
ENSG00000233653.3   0   0   0   0   0   0   0
ENSG00000250575.1   1   11  17  5   5   8   17
ENSG00000278757.1   0   0   2   0   1   0   7
ENSG00000230021.9   0   25  30  27  9   8   11
ENSG00000235146.2   0   0   0   0   0   0   0
ENSG00000225972.1   0   7   3016    11  5   14  5
ENSG00000225630.1   1   113 194 76  91  47  94
ENSG00000237973.1   25  1037    9767    1160    886 1321    5220

```

What I would like to do is only keep the ~top 50 highly variable genes:

I thought it would be like this: `dif_genes <- order(apply(df, 1, var), decreasing=TRUE)[1:50]` but that leaves me with a vector containing variances without genes.

# Answer
> 3 votes
One way is adding a new `var` column to your original data frame and sorting the updated data frame by the newly introduced column:

```
df$var <- apply(df, 1, var)
df[order(df$var, decreasing=TRUE),]

```

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


title: How to improve this WGCNA analysis

text:
# Question
Title: How to improve this WGCNA analysis
This is my first time working on using WGCNA on a microarray dataset. One of the major problems I am facing is merging close modules which is not really working well. I have quantile normalised the data before working on it. I have put my code below and link to some of the plots generated during this process. I was wondering can somebody advice on what could be improved upon in this analysis? Link to plots: !\[plots\] https://i.stack.imgur.com/dF6NA.jpg I have a total of 767 samples and 13466 probes.

```
powers = c(c(1:10))
## Scale free topology since it can find hubs
sft = pickSoftThreshold(bryois_norm_keep_use, powerVector = powers, verbose = 5)
par(mfrow = c(1,2));
cex1 = 0.9;
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
labels=powers,cex=cex1,col="red");
abline(h=0.90,col="red")
# Mean connectivity as a function of the soft-thresholding power
plot(sft$fitIndices[,1], sft$fitIndices[,5],
xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
# Turn data expression into topological overlap matrix
power=sft$powerEstimate
# Even though estimate suggests 3, trying ahigh number
# Turn adjacency into topological overlap matrix (TOM)
adjacency <- adjacency(bryois_norm_keep_use, power = 6)
TOMadj <- TOMsimilarity(adjacency,verbose = 5)
dissTOMadj <- 1- TOMadj
# Clustering using TOM
# Call the hierarchical clustering function
hclustGeneTree <- hclust(as.dist(dissTOMadj), method = "average")
# Plot the resulting clustering tree (dendogram)
sizeGrWindow(12, 9)
plot(hclustGeneTree, xlab = "", sub = "",
main = "Gene Clustering on TOM-based disssimilarity",
labels = FALSE, hang = 0.04)
# Make the modules larger, so set the minimum higher
minModuleSize <- 30
# Module ID using dynamic tree cut
dynamicMods <- cutreeDynamic(dendro = hclustGeneTree,
distM = dissTOMadj,
deepSplit = 2, pamRespectsDendro = FALSE,
minClusterSize = minModuleSize,verbose = 5)
table(dynamicMods)
# Convert numeric lables into colors
dynamicColors <- labels2colors(dynamicMods)
table(dynamicColors)
# Plot the dendrogram and colors underneath
sizeGrWindow(8,6)
plotDendroAndColors(hclustGeneTree, dynamicColors, "Dynamic Tree Cut",
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05,
main = "Gene dendrogram and module colors")
# Calculate eigengenes
dynamic_MEList <- moduleEigengenes(bryois_norm_keep_use, colors = dynamicColors)
dynamic_MEs <- dynamic_MEList$eigengenes
library(dynamicTreeCut)
# Calculate dissimilarity of module eigengenes
dynamic_MEDiss <- 1-cor(dynamic_MEs)
dynamic_METree <- hclust(as.dist(dynamic_MEDiss))
# Plot the hclust
sizeGrWindow(7,6)
plot(dynamic_METree, main = "Dynamic Clustering of module eigengenes",
xlab = "", sub = "")
######################## MERGE SIMILAR MODULES
dynamic_MEDissThres <- 0.95
# Plot the cut line
#abline(h = dynamic_MEDissThres, col = "red")
# Call an automatic merging function
merge_dynamic_MEDs <- mergeCloseModules(bryois_norm_keep_use, dynamicColors, cutHeight = dynamic_MEDissThres, verbose = 5)
# The Merged Colors
dynamic_mergedColors <- merge_dynamic_MEDs$colors
# Eigen genes of the new merged modules
mergedMEs <- merge_dynamic_MEDs$newMEs
mergedMEs
table(dynamic_mergedColors)
sizeGrWindow(12,9)
plotDendroAndColors(hclustGeneTree, cbind(dynamicColors, dynamic_mergedColors),
c("Dynamic Tree Cut", "Merged dynamic"),
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05)

```

# Answer
> 3 votes
I took a quick look and don't see anything wrong with the module merging. You could lower the merging threshold a bit - 0.95 is very high, I would use 0.9 or 0.85. Other than that though, really nothing seems out of the ordinary.

---
Tags: wgcna
---


# Question
Title: Make mpileup file with several BAM files
I'm new to this business, so my question may seem as dumb. I want to make an mpileup file using a reference sequence and several (two to begin with) aligned sorted BAM files. And I want the output file to have certain easy-to-read tab-delimited formatting (which is helpful for subsequent analyses of the output file):

field 1 - chromosome,
field 2 - position,
field 3 - reference base,
field 4 - base(s) in SAMPLE1.BAM for position in field 2,
field 5 - number of reads supporting bases in field 4,
field 6 - base(s) in SAMPLE2.BAM for position in field 2,
field 7 - number of reads supporting bases in field 6.

Something like shown below:

```
field1   field2   field3   field4   field5   field6   field7
Chr1     15967    T        T,A      15,17    T,_      25,_

```

where _ means no alternative (or second) base (or zero reads) for this position in the BAM file. '0' or empty space (or nothing) instead of _ is ok. It's ok if information in fields 4-7 is shown not in a single field but in several fields.

Which command is suitable for such a purpose? Thank you.

# Answer
> 3 votes
Have you tried using `samtools mpileup in.bam [in2.bam [...]]`? Here's information about the output format from the manual page:

> *The first three columns give the position and reference:*
>
> * Chromosome name.
> * 1-based position on the chromosome.
> * Reference base at this position (this will be “N” on all lines if -f/--fasta-ref has not been used).
>
> *The remaining columns show the pileup data, and are repeated for each input BAM file specified:*
>
> * Number of reads covering this position.
> * Read bases. This encodes information on matches, mismatches, indels, strand, mapping quality, and starts and ends of reads.

This format seems nearly identical to what you've asked for, so hopefully is close enough for you to tweak your downstream analyses.

---
Tags: sam, mpileup
---