title: GROMACS RMSF Result Analysis Problem

text:
# Question
Title: GROMACS RMSF Result Analysis Problem
This is RMSF for protein from GROMACS after 100ns simulation-

I have few questions-

1. Residue numbers are starting from 1650 to 1850, has it ignore LIG and HOH as residues?
2. There are three lines, so are there three chains and their residue number is same?
3. Are the straight lines connecting these chains or this is kind of loop something?
4. Since I am seeing 2nm or 20A fluctuation, is the protein quite disturbed?

# Answer
> 2 votes
1. The residue numbering in a PDB file mainly depend on the person who deposited it. There are usually reasons for how things are numbered, but ultimately you should look at PDB file, and check the details on how the structure was obtained. To give you an example, during my PhD I worked on a protein that had some missing residues at the N-terminus. That was because the N-term was very unstable and they didn't manage to crystallise it properly.
2. You need to check the PDB, but it is possible that the 3 chains have the same residue numbering (A1-200, B1-200). This is common in structures where you have multiple identical sub-units (e.g. dimers).
3. The lines are merely an artefact. Those lines are the attempt of the engine to draw a single line for all the chains. You go from residue 1850 of chain 1 to residue 1650 of chain 2, and so on. Some tools handle this behaviour, others don't. I would suggest you to parse the PDB yourself and plot (e.g. `ggplot`) each chain as a layer of the same figures .
4. It is indeed, but only you can tell if this is good or bad. It really depends on how you calculate the RMSF. Was it on the backbone or on the side-chains? How is your dynamics? High temperature? Does you protein partially unfold? are the chains separating?

---
Tags: protein-structure, gromacs, molecular-dynamics
---


title: WGCNA hub genes error

text:
# Question
Title: WGCNA hub genes error
I am trying to find hub genes for modules made using WGCNA:

```
hub = chooseTopHubInEachModule(datExpr, moduleColors)

```

But I am getting error:

```
Error in hubs[m] <- colnames(adj)[hub] : replacement has length zero

```

Thank you for the help1

# Answer
> 1 votes
The most likely culprit here is that your expression data (the datExpr input to chooseTopHubInEachModule) has no colnames. Check that your data do have appropriate column names. Also check that the module color vector you supply has at least actual module (type table(moduleColors) and check that you have at least one module apart from "grey" or 0).

---
Tags: wgcna
---


title: dbSNP GRch37 Download

text:
# Question
Title: dbSNP GRch37 Download
I am still new to this.

I am trying to download the GRch37 dbSNP data from the NCBI dbSNP Browser using the "Send to:" button. But whenever I select

"Send to:" \> "File" \> "Tab" \> "Default order"

The coordinates of my positions will be GRch38 based. Does anyone knows how can I download the dbSNP data with GRch37 coordinates?

Thanks

Ken

# Answer
> 2 votes
I see no option to specify the reference version using the NCBI dbSNP browser, it is set to download the latest built information. I think this is not the standard tool to download dbSNP data across many sites, but rather to browse it and visualize it for specific variation.

If you are only interested in getting info for a smaller range of positions, I recommend you to use the UCSC browser, go to Table Browser, -\> assembly GRch37, -\> group variation, -\> track all SNPs, -\> table snp51 -\> have a look at "describe table schema" to see what will you get, select your genomic range, and click on get output.

Alternatively and if you have a lot of space available in a server, simply go to the NCBI dbSNP FTP site, and download the latest GRch37 dataset (GCF\_000001405.25) with all the entries. Make sure you read the README for further info about the files and data.

---
Tags: public-databases, phylogenetics
---


title: How to Read SCF file in Python?

text:
# Question
Title: How to Read SCF file in Python?
Is there any way that I can read SCF file in python like in R using sangerseqR, I have tried with Biopython, it seems it does not support this format.

# Answer
> 1 votes
It was done using sangerseqR with rpy2 in python

```
    import rpy2.robjects as r
    from rpy2.robjects.packages import importr
    utils = importr('utils')
    utils.install_packages('sangerseqR', repos="https://git.bioconductor.org/packages/sangerseqR")
    utils.chooseBioCmirror(ind=1) # select the first mirror in the list
    utils.install_packages('sangerseqR')

    sangerseqR = importr('sangerseqR')
    Trace = sangerseqR.readsangerseq( file="1I1_F_P1815443_047.scf") #scf sequence file is now loaded

```

sangerseqR Documentation(pdf)

---
Tags: python, ngs, sequence-analysis, data-preprocessing
---


# Question
Title: Mapping unique GO term description given a specific GO ids
*This question was also asked on Biostars*

I have a list of GO ids and I want to find out a unique term description such that if I provide say 200 GO IDs I will give 200 specific GO terms. The code snippet I am using is given below. Although it kind of works, for one given GO ID I get several terms like for my 200 GO ids I get 2577 terms whereas I am interested in a specific GO term. The output that I want is like this GO:0007254 = JNK cascade. I will like to mention that I have also cross-posted the query in biostars and got no reply for the last 2 days and hence I am posting the query on Stack exchange bioinformatics. Any help will be useful.

```
library(biomaRT)
mart <- useMart(biomart = "ENSEMBL_MART_ENSEMBL", dataset = "drerio_gene_ensembl", host = "asia.ensembl.org")
v1 <-c(target)   # The list of GO ids that I want to pass as a list.
go_list <- getBM(attributes=c("go_id", "name_1006", "namespace_1003"),
              filters = "go",
              values = c("v1"),
              mart=mart, uniqueRows = TRUE)

```

# Answer
```
library(biomaRt)
mart <- biomaRt:::useMart(biomart = "ENSEMBL_MART_ENSEMBL", dataset = "drerio_gene_ensembl", host = "asia.ensembl.org")

v1 <- c("GO::0009055","GO:0016020","GO:0016491","GO:0006122","GO:0008121")

go_list <- biomaRt:::getBM(attributes=c("go_id", "name_1006", "namespace_1003"),
              filters = "go",
              values = v1,
              mart=mart, uniqueRows = TRUE)

```

so it is likely the "v1" that should not be in quotes

```
head(go_list)
        go_id                                                   name_1006
1 GO:0016020                                                    membrane
2 GO:0005739                                               mitochondrion
3 GO:0005743                                mitochondrial inner membrane
4 GO:0070469                                                 respirasome
5 GO:0006122 mitochondrial electron transport, ubiquinol to cytochrome c
6 GO:0001525                                                angiogenesis
      namespace_1003
1 cellular_component
2 cellular_component
3 cellular_component
4 cellular_component
5 biological_process
6 biological_process

```
> 2 votes

---
Tags: genomics
---


# Question
Title: How does cromwell and AWS batch decide which compute environment (on-demand or spot) to use?
For processing genomics workflows I make use of widdle and cromwell with an AWS back-end. I have followed all examples on the AWS cromwell tutorial. In my AWS batch console I can see there are 2 job queues (priority-cromwell and default-cromwell) and both queues can submit jobs to 2 compute environments (ondemand-cromwell and spot-cromwell). I have however no idea how I can control which jobs are put in which queue and are executed in which compute environment. I'd like to make optimal use of the spot-instances, how can I do that? Do I have to add something in the runtime block of my widdle task definitions, should I include something in the cromwell configuration file, or are for example only scattered jobs executed on spot instances?

# Answer
> 1 votes
Cromwell will send jobs to the AWS Batch queue ARN defined in the `cromwell.conf` file in the Cromwell server. To use Spot instances you can create an AWS Batch Compute environment configured to use Spot instances and create a queue for that compute environment. You can then set the arn of that queue in the configuration file or override the configured value using the standard mechanisms of overriding Cromwell’s configuration (workflow options, runtime parameters etc). If you change the configuration file you will need to restart the server.

Batch queues can point to more than one compute environment. Each will have a priority. Jobs in a queue will be sent to the highest priority environment until that environment has reached its max vCPU capacity, at that point jobs will be sent to the next environment.

---
Tags: workflow-management, aws
---


# Question
Title: Same region/region but different annotation
I have done the differential accessibility analysis then annotated the region with chipseeker.

my output for one of the many region and question regarding this

```
Peak log2FoldChange         padj seqnames    start      end width strand        annotation geneChr geneStart
105  Peak32946      -1.992474 9.197721e-19    chr11 14078815 14079558   744      * Distal Intergenic      11  13962689
5445 Peak32938      -2.473959 4.878995e-04    chr11 13961467 13961818   352      *          Promoter      11  13962689
      geneEnd geneLength geneStrand geneId      transcriptId distanceToTSS         ENSEMBL SYMBOL  GENENAME
105  14267884     305196          1  10418 ENST00000576479.3        116126 ENSG00000262655  SPON1 spondin 1
5445 14267884     305196          1  10418 ENST00000576479.3          -871 ENSG00000262655  SPON1 spondin 1

```

So the same region is marked as Distal Intergenic and Promoter. Next I have the same comparison RNA seq expression data. So now my objective is to find those genes+region which are expressed as well as deferentially accessible.

So for example here in this case in RNA seq data there is only one `SPON1 gene` so if I have to show it associated region how to put that,shall I consider promoter or the distal intergenic or consider which is more significant.

# Answer
You have two different peaks at different positions with different annotations. This makes complete sense, you're just misreading it.
> 2 votes

---
Tags: r, atac-seq
---


# Question
Title: How can I calculate the distances between two specific residues of a protein from a PDB file?
Suppose, we have a protein (`16PK`).

We are considering a 5-residue segment/window.

If we assume `ASN10` to be the `i`th residue, how can I calculate the following distances in **Python**?

1. from \[i+2\] to \[i-2\]
2. from \[i+1\] to \[i-2\]
3. from \[i+2\] to \[i-1\]

Next major question is, what would be the distances like if my segment/window is 7 or 9 residues long?

# Answer
Your $\[i,j,k\]$ coordinates in the PDB file are coordinates. You can readily compute Euclidean distances between any pair of coordinates as described in the wiki page:

$d(p,q)=\sqrt{(p\_1-q\_1)^2 + (p\_2-q\_2)^2 + (p\_3-q\_3)^2}. $

This resource suggests one such approach (code reproduced below in case that site stops existing):

```
import Bio.PDB
import numpy

def calc_residue_dist(residue_one, residue_two) :
    """Returns the C-alpha distance between two residues"""
    diff_vector  = residue_one["CA"].coord - residue_two["CA"].coord
    return numpy.sqrt(numpy.sum(diff_vector * diff_vector))

```
> 2 votes

---
Tags: python, proteins, pdb
---


# Question
Title: Counting the co-occurrence of variants A and B in an aligned sequencing read
I need some help getting started on this project.

To simplify we want to be able to quantify the occurrence of 3 variants on each sequencing read in an alignment file as a proxy measurement for recombination.

The possible combinations we have are A, B or AB.

For example, variant A is at position 100 and variant B is at position 300 of the reference and AB would have both.

We have alignment files against a reference sequence, and then I get stuck. I could technically look by eye in tablet or IGV for example but is there a way to extract each read then query by expected base at the reference position then bin these reads into the categories; A, B, AB and unclassified then calculate the total occurrences of each?

Thanks for your time and suggestions!

# Answer
Related question: Access base aligned to particular reference position

After an initial effort posted here, I ended up rewriting this and testing it a little, and posted the script here. It still has some drawbacks (handling indels intelligently) but now it more or less works.
> 1 votes

---
Tags: sam, samtools, mutations, variants
---


# Question
Title: How to convert vcf to tbed/tfam and again convert it back to tbed/tfam to vcf using plink
When converting vcf to tbed/tfam format, I am getting error like

```
unrecognized flag<'--transpose'>

```

Code I used:

```
plink2 --vcf sample_vcf_file.vcf --recode --transpose --out outputfile

```

# Answer
> 1 votes
Correct command is

```
plink2 --vcf sample_vcf_file.vcf --recode transpose --out outputfile

```

vcftools also can convert vcf into tped/tfam

---
Tags: plink
---


# Question
Title: Translate all reads in a .fastq into protein sequence from deep mutational scanning experiment
I'm working with paired-end NGS reads from an Illumina platform. The sequences I have are all of the same gene, but have one or more substitution mutations each.

Here is a rough workflow for generating the reads:

1. Expression of a library of sequence variants for a heterogeneous protein of interest in phage.
2. Selection of high-fitness variants according to a custom assay.
3. Extraction of the variant pool via amplification of the variable region of the gene of interest.
4. Preparation of Illumina sequencing library from the amplicon pool.
5. Deep sequencing of the amplicon pool.

I've already combined the paired-end reads with USearch and manually removed sequences with subpar quality scores, but I now need to translate all of the remaining DNA sequences into amino acid sequences. I'm uncertain of where to begin addressing this issue, and I was hoping someone could point me in the direction of an application or python package that could assist me.

My goal: **Count the abundances of each sequence variant within the Illumina data to identify which of the *protein* sequences are the most fit.**

# Answer
Based on back and forth in the comments, I think that you might be interested in a deep mutational scanning workflow, for example `dms_tools2` or `enrich2`.

These are built specifically for the kind of experiment that you are interested in, if I understand correctly. See attached figure from dms\_tools2:
> 2 votes

---
Tags: python, proteins, phylogenetics, fastq
---


# Question
Title: What does the following diagram represent?
Suppose, the following is a diagram of a protein's polypeptide chain:

What does this diagram represent?

1. What are the letters A, E, M, W, L, N, S, etc. represent? (I suppose these are amino acid codes. Right?)
2. Which one is the c-alpha atom here?
3. What are the sequence numbers `i-2`, `i-1`, `i`, `i+1`, `i+2`, etc. represent? Why there are negative and positive increments from `i`?

The following is a related diagram:

# Answer
> 2 votes
1. Yes, the letters are amino acids
2. It's a schematic with a whole residue as a circle —possibly from Voet and Voet.
3. Residue `i` is the one in focus, the reference. `i-1` is the preceding residue. Have a gander at en.wikipedia.org/wiki/Turn\_(biochemistry) where this type of relative notation is important.

The diagram shows how secondary structure (helices, sheets etc.) affect the distance between the residues. In a helix a residue 3 positions down will be close, but not in a sheet.

However, for assignment of secondary structure it may be easier to look at the φ and ψ angles:

* α helix: phi=-57.8, psi=-47.0
* π helix: phi=-57.1, psi=-69.7
* 3.10 helix: phi=-74.0, psi=-4.0
* β sheet: phi=-139, psi=+135

Although these have a range and you can get π bulges within α helices —important for conformation switching for pore opening (capsaicin binding to TRPV1 is a classic example).

---
Tags: proteins, protein-structure
---


# Question
Title: no expression for hub genes identified in wgcna
I tried to identify intramodular hub genes in `WGCNA` analysis. I used:

```
hub = chooseTopHubInEachModule(datExpr, moduleColors)

```

But when I checked the hub genes for the modules, I found that the hub genes are expressed only in less than 5 out of 300 samples. I did not filter the gene expression data for any low expressed genes in certain number of samples. How does this make biological relevance if the hub genes are not expressed in most of the samples?

Thank you!

# Answer
> 1 votes
It could be that your modules group together genes (over)expressed in only certain small groups of samples that may be biologically relevant. You would need to check whether the samples in which the hubs (and presumably the modules) are up share some sort of a characteristic.

In general though, there are good reasons why authors of WGCNA recommend filtering out non- and low-expressed genes - correlations tend to not be a good co-expression measure when counts are very low, and technical effects can be much more pronounced. It is also recommended to spend a bit of time on preprocessing to make sure the data set does not have strong global drivers that are technical or uninteresting.

---
Tags: wgcna
---


# Question
Title: How do I tell that a MinION run has finished based on the files produced?
I want to automate some scripts which should run after a MinION run completes. Is there a file which is produced by MinKNOW which will enable me to tell that the run has finished?

# Answer
You can just look for the presence of a file called "final\_summary\_\*.txt" in the run folder. From Matt Loose.
> 1 votes

---
Tags: nanopore, minion
---


# Question
Title: How to get a FASTA format file with the DNA sequences of all annotated genes?
I am analysing Pyrococcus Furiosus DNA sequencing data by considering data published here in NCBI. When I click on "Send to"\>"Gene Features"\>"FASTA format" I download a file that has the sequences of genes of this organism but I realized that this file has some sequences that are doubled (that is the sequence is present twice with different identifiers)... is there a way to get a file with all genes annotated and their respective DNA sequences without double sequences and with all sequences in the 5' to 3' direction without reverse complement in a well "ordered" way ? In NCBI (in the link I reported before) it indicates 2,128 genes so I would like a file with all these 2,128 genes annotated and their respective DNA sequences in FASTA format. Do you know if there is an other website or an other place in NCBI in which I can search to get this kind of file?

EDIT: if you follow the path "Send to"\>"Gene Features"\>"FASTA format", now it seems that there are just 2,128 genes so this problem is vanished in this case. Anyway, I would appreciate the answer for the future and for the reverse complement gene sequences that are still there.

Thank you very much.

# Answer
From this SO post. You can use seqkit to remove duplicate sequences with the command below:

```
seqkit rmdup -s < sequences.txt > out.fa

```

> The `rmdup` option removes duplicates, and the `-s` option calls duplicates on the basis of sequence, ignoring differences in headers.
> 1 votes

---
Tags: fasta, sequence-analysis, public-databases, phylogenetics, public-dataset
---


# Question
Title: Accessing expression data in an ExpressionSet
I used `rma` from the oligo package to normalize CEL files, returning an object of class `ExpressionSet`. How do I access the actual expression data?

```
> data = read.celfiles(names)
> data2 = oligo::rma(data)
class(data2)
[1] "ExpressionSet"
attr(,"package")
[1] "Biobase"
> data3 = as.matrix(data2)
Error in as.vector(data) :
  no method for coercing this S4 class to a vector

```

# Answer
The class of your object which the `rma` function returns is an `ExpressionSet` which is a standardized data structure in Bioconductor/R from the Biobase package.

It stores both expression (a numeric matrix) data as well as column- and row annotation data. You have to use the getter function `exprs()` to access the expression data slot, so `exprs(data2)` which is then of class `matrix`, link to documentation.
> 2 votes

---
Tags: bioconductor, microarray
---


# Question
Title: Transfer annotations from one genome assembly to another
I've got an annotated draft genome assembly made using short and long read strategies. I've also done optical mapping to stitch some of the contigs together, with the goal being chromsome-length assemblies. How do I annotate the output .fasta from the optical mapping assembly using the .gff of the draft assembly that was used as input? I also have a .agp file that describes what sequence fragments ended up where in the output.

# Answer
> 2 votes
You should use a liftover approach. You could either do this based on exact changes, which might be a bit tedious given the presumably large number of join operations you've done (especially if you've e.g. broken contigs).

Therefore it might be more convenient to do it on the basis of homology.

liftoff is one simple tool for doing this with a gff file and your 2 fasta files. It seems to mostly work ok.

---
Tags: fasta, sequence-annotation, genomics
---


# Question
Title: Create BED file from list with gene symbols
Is there a way to generate a BED file of a list of gene symbols for a given genome?

# Answer
> 4 votes
Via Gencode and BEDOPS `convert2bed`:

```
$ wget -qO- ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.annotation.gff3.gz \
    | gunzip --stdout - \
    | awk '$3 == "gene"' - \
    | convert2bed -i gff - \
    > genes.bed

```

Replace the FTP path with that needed for your genome of interest. Feel free to edit your question if you need help with that, by specifying your given genome.

You can modify the `awk` statement to get exons, by replacing `gene` with `exon`.

BEDOPS: https://github.com/bedops/bedops

Once you have the genes in BED format, you can `grep` the file for them, using a text file of matches:

```
$ grep -Fwf genes_of_interest.txt genes.bed > genes_of_interest.bed

```

# Answer
> 3 votes
You could also get your genes from UCSC. Here's how:

---
Tags: fasta, genome, bed
---


# Question
Title: Why are my Nextflow processes not executing in parallel?
I have written a Nextflow script with three process:

* The first process takes a pair of `fastq` files and aligns with reference genome. The process writes the resulting `SAM` file into `sam channel`.
* Second process takes input from the `sam channel` and creates a `BAM` file from it, and writes it into `bam channel`.
* Final process reads from `bam channel` and sorts the `BAM` files.

Following is the entire code:

```
#!~/bin nextflow

reads_ch = Channel.fromFilePairs('raw/*_{1,2}.fastq.gz', flat: true)
ref_genome = file('reference/genome/index/ref_genome')

process mapToRef {
    memory '20 GB'
    cpus 16
    input:
    tuple val(sample_id), file(read1), file(read2) from reads_ch

    output:
    tuple val(sample_id), file("${sample_id}.sam") into sam_ch

    script:
    """
    ~/biotools/bwa/bwa mem -t 60 ${ref_genome} ${read1} ${read2} > ${sample_id}.sam
    """
}

process samToBam {
    memory '20 GB'
    cpus 16
    input:
    tuple val(sample_id), file("${sample_id}.sam") from sam_ch

    output:
    tuple val(sample_id), file("${sample_id}.bam") into bam_ch

    script:
    """
    ~/biotools/SAMTOOLS/samtools-1.12/samtools view -S -b ${sample_id}.sam > ${sample_id}.bam
    """
}

process sortBam {
    memory '20 GB'
    cpus 16
    input:
    tuple val(sample_id), file("${sample_id}.bam") from bam_ch

    output:
    tuple val(sample_id), file("${sample_id}.sorted.bam") into sorted_bam_ch

    script:
    """
    ~/biotools/SAMTOOLS/samtools-1.12/samtools sort ${sample_id}.bam -o ${sample_id}.sorted.bam
    """
}

```

All three processes successfully execute their commands. However, they are executed sequentially. For example, the first process `mapToRef` runs **sequentially** for all read pairs, followed by `samToBam`, which too, runs sequentially after all `mapToRef`s have finished execution. I was under impression that parallelization would "work" out of the box, but I must be doing something wrong. I am a complete beginner in Nextflow and would greatly appreciate any help regarding this.

**EDIT**

Output of `reads_ch.view()`:

```
[ERR2512393, /home/work/raw/ERR2512393_1.fastq.gz, /home/work/raw/ERR2512393_2.fastq.gz]
[ERR2512394, /home/work/raw/ERR2512394_1.fastq.gz, /home/work/raw/ERR2512394_2.fastq.gz]
[ERR2512391, /home/work/raw/ERR2512391_1.fastq.gz, /home/work/raw/ERR2512391_2.fastq.gz]
[ERR2512392, /home/work/raw/ERR2512392_1.fastq.gz, /home/work/raw/ERR2512392_2.fastq.gz]
[ERR2512390, /home/work/raw/ERR2512390_1.fastq.gz, /home/work/raw/ERR2512390_2.fastq.gz]

```

# Answer
> 4 votes
I have tested it on my end with the following script:

```
#!/usr/bin/env nextflow

reads_ch = Channel.fromFilePairs('/root/test_nextflow/*_{1,2}.fastq.gz', flat: true)

process mapToRef {
  memory '1 GB'
  echo true
  cpus 1

  input:
  tuple val(sample_id), file(r1), file(r2) from reads_ch

  output:
  tuple val(sample_id), file("${sample_id}.sam") into sam_ch

  script:
  """
  echo bwa mem ${r1} {r2}
  if [ ${sample_id} == "ERR2512393" ]; then sleep 10; fi
  touch ${sample_id}.sam
  """
}

process samToBam {
  memory '1 GB'
  cpus 1
  echo true

  input:
  tuple val(sample_id), file(sam) from sam_ch

  output:
  tuple val(sample_id), file("${sample_id}.bam") into bam_ch

  script:
  """
  echo samtools view -S -b ${sample_id}.sam
  touch ${sample_id}.bam
  """
}

process sortBam {
  memory '1 GB'
  cpus 1
  input:
  tuple val(sample_id), file(bam) from bam_ch

  output:
  tuple val(sample_id), file("${sample_id}.sorted.bam") into sorted_bam_ch

  script:
  """
  echo samtools sort ${sample_id}.bam
  touch ${sample_id}.sorted.bam
  """
}

```

It's producing the following result:

As you can see, the two processes are running in parallel (`2 of 2`)

It could be that yours are limited by the resources? If your PC doesn't have 40GB of memory and 32 CPUs, it will run them sequentially.

**EDIT**: By adding `sleep 10` to `mapToRef`, we can clearly see both of them are running in parallel now:

**EDIT 2**: Added condition to "sleep" on one of the samples:

**Final EDIT**: be mindful of the resources you allocate to your processes. You could use `${task.cpus}` and `${task.memory.toGiga()}` to facilitate the configuration of your processes so that they are not different from what you have defined and what is actually used. See CPUs

---
Tags: nextflow
---


# Question
Title: How to see effect of point mutation in PyMOL?
I am trying to model how a specific point mutation would affect the protein structure. So far, I've figured out how to create the point mutation in PyMOL using the mutagenesis wizard. I'm not sure how to take the next step toward actually visualizing the consequences of the mutation though (i.e. would it disrupt a helix, etc). Is there a way to do this in PyMOL, or would I need other software? Thank you!

# Answer
> 3 votes
A common metric used to assess if a mutation destabilises the structure is the difference in Gibbs free energy. PyMOL does not give a ∆∆G (difference in Gibbs free energy, kcal/mol), but a very crude steric clash nor does it repack neighbouring sidechains out of the way, or alter the backbone, so unfortunately is not the way forward. However, if one must there is a way (see SIMBAI below).

## Beyond destabilisation

Destabilisation predictors have a ~60% or lower success at predicting if something is pathogenic (ref). This is because the effect a mutation has on a protein is rather complicated as you can have destabilisation, but also aggregation, loss of regulation, mislocalisation, loss of binding to another protein, loss of binding to ligand, loss of activity (catalysis, inability to change conformation) etc etc. Then you can get curious subcases, such as loss of activity that result in the sequestration of a binding partner or homologues in a complex resulting in a worse phenotype than destabilisation or truncation.

## Modelling

The behaviour of a protein can be assess by the energy of the system using forcefields, which are multi-term equations that describes the atomic interactions (attraction, repulsion etc.). The result is the Gibbs free energy mentioned above, which is a potential, so negative is good. Kinetic energy is release when the unfolded protein "rolls down" down the potential energy funnel —ProTherm dataset (an empirical dataset of ∆∆Gs) and the tool SDM give this value inverted as the negative can be confusing.

An uphill change in the Gibbs free energy (∆∆G \> 0) generally correlates with a destabilisation of a protein —aggregation and solubility being the caveat in the "generally". The destabilisation a single state (static conformation) can be assessed. `gmx energy` in Gromacs, FoldX or Rosetta can be used for this.

For altered ability to switch conformation steered MD simulations are required, but generally the ∆∆G at different states can be used.

For the non-destabilisation options, https://www.phosphosite.org/ and http://elm.eu.org/ are very useful.

I have been working on a server that gives a ∆∆G value and note if there are any nearby post-translational modifications, gnomAD or changes in potential linear motifs etc.: Venus (unpublished). However, there are many ∆∆G calculation servers out there and are compared in many papers (example), but it should be said that the dataset used in these comparisons is always ultra-biased despite how much is tried. Missense by the Phyre2 group (Michael Sternberg) give detailed information of the neighbourhood changes the mutation causes. In terms of extra info, Miscast gives Uniprot info.

## SIMBA-I via PyMOL

In terms of PyMOL, the SIMBA-I algorithm can be used by getting the RSA from it —I would not recommend it as it may underestimate the severity of unusual cases. But if one must, in PyMOL one can paste the following into the console with the wild type variant:

```
# defined by you so change accordingly
from_residue = '' # single letter
to_residue = '' # single letter
resi =
chain = ''

## data
maxASA = {'A': 121.0, 'R': 265.0, 'N': 187.0, 'D': 187.0, 'C': 148.0, 'E': 214.0, 'Q': 214.0, 'G': 97.0, 'H': 216.0,\
          'I': 195.0, 'L': 191.0, 'K': 230.0, 'M': 203.0, 'F': 228.0, 'P': 154.0, 'S': 143.0, 'T': 163.0,\
              'W': 264.0, 'Y': 255.0, 'V': 165.0}
volumes = {'A': 88.6, 'R': 173.4, 'N': 114.1, 'D': 111.1, 'C': 108.5, 'Q': 143.8, 'E': 138.4, 'G': 60.1,\
               'H': 153.2, 'I': 166.7, 'L': 166.7, 'K': 168.6, 'M': 162.9, 'F': 189.9, 'P': 112.7, 'S': 89.0,\
               'T': 116.1, 'W': 227.8, 'Y': 193.6, 'V': 140.0}
hydrophobicities = {'A': 8.1, 'R': 10.5, 'N': 11.6, 'D': 13.0, 'C': 5.5, 'Q': 10.5, 'E': 12.3, 'G': 9.0,\
                    'H': 10.4, 'I': 5.2, 'L': 4.9, 'K': 11.3, 'M': 5.7, 'F': 5.2, 'P': 8.0, 'S': 9.2, 'T': 8.6,\
                    'W': 5.4, 'Y': 6.2, 'V': 5.9}
cmd.set('dot_solvent', 1)
cmd.set('dot_density', 3)
sasa = cmd.get_area(f'resi {resi} and chain {chain}')
rsa = sasa / maxASA[from_residue]
offset = -1.6
rsa_factor = +2.2
dV_factor = +0.59
dH_factor = -0.29
rsa_dV_factor = -0.54
rsa_dH_factor = +0.49
# calculate
dV = (volumes[to_residue] - volumes[from_residue]) / 100
dH = (hydrophobicities[to_residue] - hydrophobicities[from_residue])
simba = offset + rsa_factor * rsa + dV_factor * dV + dH_factor * dH + rsa_dV_factor * rsa * dV + rsa_dH_factor * rsa * dH
print -simba  # negative stabilised, positive destabilised

```

## Human reasoning

However, the major problem with all the tools mentioned is that human reasoning is often required to read the literature and make a solid hypothesis.

A resource I cannot stress enough, for human mutations, its utility is gnomAD, which contains *phenotypically neutral* mutations: for example, one really ought to not claim a homozygous missense is pathogenic because it is causing a loss of function by destabilising the protein, yet there are homozygous truncations roaming around in the human population.

The paper describing the crystal structure used and those of the homologues may point out key features which would missed by automated servers, such as secondary-structure elements that are speculated to perform certain roles. In fact, without doing any calculations seeing in PyMOL if a mutation causes a clash in the core of a protein, alters the active site, or changes a residue to a smaller one in a region known to alter conformation upon a given condition, is a lot more informative even if not quantitative...

# Answer
> 1 votes
You could run a molecular dynamics simulation for example with Gromacs or Charm. The problem with pymol is that the structure is static and you don't know how much tension that point mutation creates. First of all you would run an energy minimization to remove bad non-physical overlaps between atoms, then add solvent molecules, and then simulate your system first with position restraints for equilibration and then without for 1ns to 100ns+.

Even then you cannot be sure that what comes out is the 'final' result as sometimes an induced fold change can take much longer. It depends on the system and the energy barriers included and the time scales that you are interested in. If the protein is large it may have many different physiologically relevant conformations and the point mutation might only change the equilibrium distribution of those states or totally prevent a transition.

In case you have more than one crystal structure (in different conformations) you may introduce these point mutations in all of them and simulate them for a while. But it is not that trivial to setup such a simulation you should contact an expert. There are many things that can go wrong. Other software you could use are Modeller and Rosetta, if I am not mistaken.

# Answer
> 1 votes
You could give SAAPpred a try. This tool uses machine learning (in the form of `random forest`) to predict the **local** effects of a single point mutation on the protein structure.

I did my PhD in the group and I used SAAPpred to find the most (predicted) damaging mutations of G6PD, to then ran extensive MD simulations to evaluate the final stability of those mutants.

Keep in mind that the tools doesn't return a structure, but just the possible/likely effects of that substitution.

---
Tags: pymol
---


# Question
Title: What is meant by transcriptional changes executed by the cell over a time period?
I read the following line in the research paper - The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells:

> During differentiation, for example, each cell makes independent fate decisions by integrating a wide array of signals from other cells and executing a complex choreography of gene-regulatory changes. Thus individual cells can execute the same sequence of transcriptional changes over highly varying time scales.

Can anyone explain in simple terms what is meant by **"same sequence of transcriptional changes"**?

How can it be the same sequence if during differentiation, different types of cells are formed?

# Answer
> 1 votes
Not all cells will perform the same sequence of transcriptional changes. The paragraph you quote simply states that differentiation isn't a step change with all cells differentiating at the exact same moment. Rather, it's a more gradual change with cells each working on their own time scale.

---
Tags: gene, single-cell, phylogenetics, rna, transcriptome-regulation
---


# Question
Title: If regulatory proteins decide the gene regulation, is it necessary for every cell to have a different look?
Transcription factors regulate the genes in different cells. That is, different cells have different sequence of genes that are transcribed (set on or off) by the RNA polymerase. This will help different types of cells to perform their own unique operation.

So, why do the cells have to look different physically, if the their function is determined by the on/off sequence?

**NOTE:** I am a neophyte in Bioinformatics. Any help on my building my conceptual knowledge is appreciated.

# Answer
There are many more contributors to the cell differentiation that the regulators. But one that you will encounter at an higher level is the chromatin status: if it is in the open status it allows transcripts to be on and being regulated (with all the further levels of regulations such as methylation, glycosylation etc). If the chromatin is off (the state is very condensed) very few transcript can be made from that region. Check this wikipedia page here
> 2 votes

---
Tags: transcription-factors
---


title: Subset GTF file for specific genes

text:
# Question
Title: Subset GTF file for specific genes
I have a GTF file (in 'test.gtf' file):

```
chr1    HAVANA  gene    10069   18412   .   +   .   gene_id "ENSG00000223972.4"; transcript_id "ENSG00000223972.4"; gene_type "pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "pseudogene"; transcript_status "KNOWN"; transcript_name "DDX11L1"; level 2; havana_gene "OTTHUMG00000000961.2";
chr2    HAVANA  gene    11869   14412   .   +   .   gene_id "ENSG0000023172.1"; transcript_id "ENSG000023172.1"; gene_type "pseudogene"; gene_status "KNOWN"; gene_name "PCNA"; transcript_type "pseudogene"; transcript_status "KNOWN"; transcript_name "PCNA"; level 2; havana_gene "OTTHUMG00000001931.2";
chr3    HAVANA  gene    15769   19412   .   +   .   gene_id "ENSG00000433972.1"; transcript_id "ENSG000023172.1"; gene_type "pseudogene"; gene_status "KNOWN"; gene_name "USP21"; transcript_type "pseudogene"; transcript_status "KNOWN"; transcript_name "USP21"; level 2; havana_gene "OTTHUMG00000002361.2";

```

I have loaded it as a pyranges object and I want to subset it so that it only contains gene from `test_list`:

```
import pyranges as pr
import pandas as pd

gtf_file = "test.gtf"
test_list = ["PCNA","USP21", "USP1"]

gtf = pr.read_gtf(gtf_file)

sub_gtf = gtf

sub_gtf.df = gtf.df[gtf.df.gene_name.isin(test_list)]

```

But the last line throws this error:

```
File "<ipython-input-24-af0b4b5c7576>", line 1, in <module>
    irongtf.df = irongtf.df[irongtf.df.gene_name.isin(test_list)]

  File "/home/niek/.local/lib/python3.8/site-packages/pyranges/pyranges.py", line 314, in __setattr__
    _setattr(self, column_name, column)

  File "/home/niek/.local/lib/python3.8/site-packages/pyranges/methods/attr.py", line 48, in _setattr
    df.insert(pos, column_name, column)

  File "/home/niek/.local/lib/python3.8/site-packages/pandas/core/frame.py", line 4413, in insert
    value = self._sanitize_column(value)

  File "/home/niek/.local/lib/python3.8/site-packages/pandas/core/frame.py", line 4504, in _sanitize_column
    com.require_length_match(value, self.index)

  File "/home/niek/.local/lib/python3.8/site-packages/pandas/core/common.py", line 531, in require_length_match
    raise ValueError(

ValueError: Length of values (2) does not match length of index (1)

```

What would be the right way to subset the gtf file so that it only contains the rows where `gene_name` contains the gene symbols from `test_list`? Or is there another tool that I could use to do this?

# Answer
> 1 votes
It turned out to be a lot easier if I used `grep` to do this:

```
grep -w -f genes.txt gencode.v19.annotation.gtf > sub_set.gtf

```

`genes.txt` contains gene symbols on new lines.

# Answer
> 3 votes
Following up on zorbax's answer, you could read in and filter the GTF file in this way, among others:

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

import gtfparse as gp

gtf_file = "test.gtf"
test_list = ["PCNA", "USP21", "USP1"]

df = gp.read_gtf(gtf_file)
subset = df[df['gene_name'].str.contains('|'.join(test_list))]
print(subset)

```

You could further cull this subset by other attributes or columns, e.g.:

```
genes_only = subset[subset["feature"] == "gene"]

```

etc.

The `gtfparse` library is available via PyPi: https://pypi.org/project/gtfparse/ or Bioconda: https://bioconda.github.io/recipes/gtfparse/README.html

# Answer
> 2 votes
If you have your gtf file like a DataFrame you can use:

```
 df[df['gene_name'].str.contains('|'.join(test_list))]

```

---
Tags: python, gtf, pyranges
---


# Question
Title: blastp option to compare with a personal database
I'm trying to compare proteins generated with prokka (from a genome) with blastp using a virulence factor database VFDB, and I want to use some parameters obtained from this paper:

`blastp option, top hit, ≥50% identity, ≥50% alignment length, E-value of < 10 - 5`

How can I specify those parameters on the command line?

```
# create DB
makeblastdb -in VFDB.fas -dbtype prot

# run blast
blastp -query my.faa -db VFDB.fas -out results.txt .....

```

any suggestion?

# Answer
You can use these options:

```
blastp -query my.faa -db VFDB.fas \
       -perc_identity 50 -outfmt 6 \
       -evalue 10e-5 -out results.txt

```

Then select the 5th column, *qlen*, and the 13th columns, *slen*, and get the percentage higher than 50.

```
awk  -F'\t' '{ if ($NF >= 50) printf "%s\t%.2f\n", $0 , ($2/$5)*100 }' result.txt

```

Check the documentation for the full `-outfmt` column names with `blast -help`
> 5 votes

---
Tags: blastp, makeblastdb
---


# Question
Title: Why does my volcano plot show genes at the extremities of the graph?
I am confused about the genes at the top of the volcano plot. There are a lot of genes at the top, what's wrong with it? If that's normal what's the interpretation of it?

I used volcano plot tools from galaxy with p-value=0.05.

# Answer
Either your graph or your data is capped at 300. Probably the p-values themselves; whatever calculated those p-values literally can't represent a number tinier than 1x10^-300
> 0 votes

---
Tags: visualization
---


# Question
Title: What assembler is appropriate for High-Fidelity PacBio reads
What assembler is appropriate for High-Fidelity PacBio reads?

For example, canu is good for high-error PacBio reads. But what algorithm to use for HiFi reads? Would it be OK to use canu without the error-correction step?

# Answer
> 3 votes
As recently tweeted by A. Phillippy, Canu v1.9 now supports HiFi reads.

Therefore you should just read their manual and look for the implementation

:)

# Answer
> 2 votes
From the group of Heng Li there is hifiasm.

# Answer
> 1 votes
The HiFi de novo assembly tool that **PacBio** has been developing is called **IPA**.

It’s available on pbbioconda for command line use, and as a push-button pipeline in their SMRT Link software.

I looked at a few examples of runtimes of **IPA** for different genome sizes:

3 Gb human genome – ~600 CPU hours 800 Mb complex plant – ~430 CPU hours 400 Mb plant - ~100 CPU hours

Hifiasm seems to be more efficient as far as I can tell. In their paper they have a few nice tables comparing it to HiCanu, Peregrine and Falcon.

---
Tags: assembly, software-recommendation, long-reads, pacbio
---


# Question
Title: How do you set the coverage in PacBio's Sequel II?
I am reading the Whole Genome Sequencing for *de novo* Assembly Best Practices

> * Use the Sequel II or IIe System and SMRT® Cell 8M to sequence to desired coverage depth for complexity of genome
> * 10- to 15-fold coverage per haplotype recommended

First of all Wikipedia shows different definitions for the term haplotype. Which one applies here?

Second while the recommendations state to use `10- to 15-fold` coverage per haplotype I find myself confused about how I would know the haplotype beforehand.

From what I understand of haplotype-estimation I would need the sequenced data first.

How would I go about this when wanting to do a *de novo* assembly of a species that has no published assembly of any close relatives? That means the only information I would have at that point is the C-value of the species to be sequenced.

Or better, how does the Sequel II figure out the number of haplotypes for achieving the desired fold-coverage?

I hope I'm in the right SE.

# Answer
> 2 votes
They just want to know the ploidy of the sample. In a diploid sample, e.g. humans, there are two haplotypes. This is mostly important for haplotype resolution (i.e. get the two full genomes in a diploid).

Note that elsewhere in the brochure that you link it writes:

> Fast and efficient assembly of even the largest genomes with haplotype resolution of complex polyploids

If you are working with a diploid species, that translates to 20-30X coverage. That seems low to me for the older tech of CLR/Sequel sequencing, but they are the experts.

Different recommendations posted here by a service provider suggest 40-50X for structural variant analysis, which seems more reasonable.

# Answer
> 1 votes
Unless you want to reconstruct all haplotypes, you can sequence to ~30X regardless of ploidy.

---
Tags: assembly, coverage, pacbio
---


# Question
Title: Does the Protein Data Bank contain the estimated distances obtained using NMR spectroscopy?
From the Wikipedia entry on the Protein Data Bank (PDB) (emphasis mine):

> Most structures are determined by X-ray diffraction, but about 10% of structures are determined by protein NMR. When using X-ray diffraction, approximations of the coordinates of the atoms of the protein are obtained, **whereas using NMR, the distance between pairs of atoms of the protein is estimated**. The final conformation of the protein is obtained from NMR by solving a distance geometry problem.

I am currently working on an algorithm to solve distance geometry problems in the context of molecular conformations. I would prefer to test this algorithm on real data, so does the PDB contain information about the distances between pairs of atoms that were estimated from NMR spectroscopy measurements?

I checked the PDB file for the 2KB7 protein, but I only found the xyz coordinates of the atoms obtained *after* solving the distance geometry problem. Is it possible to get access to the estimated distances that were used to generate these conformations?

# Answer
> 2 votes
The data can be found here:

https://bmrb.io/search/instant.php?term=2KB7

There is no link directly in the PBD as it is a one-to-many relationship. Clicking on the first value (a shifts dataset) you get somewhere in the middle a link to the actual dataset, which is a CIF like file with the following table:

```
                              vvvvvvv
    1 . 1 1  2  2 MET HA H  1   4.070 . . 1 . . . A  2 MET HA . 18256 1
    2 . 1 1  2  2 MET C  C 13 176.670 . . 1 . . . A  2 MET C  . 18256 1
                              ^^^^^^^

```

Each row has the spectral shift in ppm (chevrons added) of an *assigned* hydrogen or carbon atom. There are also anisotropic chemical shifts and dipolar couplings datasets etc. etc. Each with a particular piece of data...

Chemical structures are "elucidated" (=solved) in different ways depending on the type of experiment.

Say you wanted only the distanced calculated by Nuclear Overhauser effect (r\_NOE), you would need to find a NOE NMR experiment.

---
Tags: protein-structure, pdb
---


# Question
Title: How to extract all sequences mapped to a transcript from Kallisto output
*This question has also been asked on Biostars*

I ran Kallisto with the --pseudobam option. How do I extract all the short reads that are mapped to a single transcript (e.g. ENST00000367969.8)? As a person without any previous SAM/BAM experience, I tried the following things without success. Really appreciate some help here.

I first tried:

```
samtools view -S pseudoalignments.bam ENST00000367969.8

```

but got the following error:

```
[E::idx_find_and_load] Could not retrieve index file for 'pseudoalignments.bam'
[main_samview] random alignment retrieval only works for indexed BAM or CRAM files.

```

I tried to index the file using:

```
samtools index pseudoalignments.bam

```

but get the following error:

```
[E::hts_idx_push] NO_COOR reads not in a single block at the end 79729 -1
[E::sam_index] Read 'A00742:184:HLTF2DSXY:4:1101:4417:1000' with ref_name='ENST00000556431.1', ref_length=4142, flags=83, pos=2098 cannot be indexed
samtools index: failed to create index for "pseudoalignments.bam"

```

# Answer
> 1 votes
Sort the BAM file again with `samtools sort` and then run `samtools index` on the result. That should fix the issue.

# Answer
> 1 votes
What Devon said. Alternatively, you could also just convert the bam to sam and then look for lines where the 3rd field (the reference the sequence is aligned to) is what you want:

```
samtools view pseudoalignments.bam | awk '$3=="ENST00000367969.8"'

```

But really, it's far better to sort and index the file as Devon suggested so you can use dedicated tools in future.

---
Tags: sam, samtools, kallisto
---


# Question
Title: Help to understand the code for dipeptide composition calculation (in python)
Dipeptide composition of a protein sequence is the number of times a particular dipeptide (e.g. Arginine-Histidine) occurs in a sequence divided by the total number of dipeptides in the sequence (which is the length of the sequence - 1)

I have found the following code for this:

```
import re

def DPC(fastas):
    AA = 'ACDEFGHIKLMNPQRSTVWY'
    encodings = []
    diPeptides = [aa1 + aa2 for aa1 in AA for aa2 in AA]
    header = ['#'] + diPeptides
    encodings.append(header)

    AADict = {}
    for i in range(len(AA)):
        AADict[AA[i]] = i

    for i in fastas:
        name, sequence = i[0], re.sub('-', '', i[1])
        code = [name]
        tmpCode = [0] * 400
        for j in range(len(sequence) - 2 + 1):
            tmpCode[AADict[sequence[j]] * 20 + AADict[sequence[j+1]]] = tmpCode[AADict[sequence[j]] * 20 + AADict[sequence[j+1]]] +1
        if sum(tmpCode) != 0:
            tmpCode = [i/sum(tmpCode) for i in tmpCode]
        code = code + tmpCode
        encodings.append(code)
    return encodings

```

There is another file (main) from where the arguments `fastas` come from which is basically a list of strings where every string is a protein sequence. I am finding it really difficult to understand the nested loops that were used:

```
    for i in fastas:
        name, sequence = i[0], re.sub('-', '', i[1])
        code = [name]
        tmpCode = [0] * 400
        for j in range(len(sequence) - 2 + 1):
            tmpCode[AADict[sequence[j]] * 20 + AADict[sequence[j+1]]] = tmpCode[AADict[sequence[j]] * 20 + AADict[sequence[j+1]]] +1
        if sum(tmpCode) != 0:
            tmpCode = [i/sum(tmpCode) for i in tmpCode]
        code = code + tmpCode
        encodings.append(code)

```

I have spent 4-5 hours on this but I am still finding it really difficult to follow. I would really appreciate if anyone can explain the steps involved in these loops. Thank you!

**Edit**: The github link to the files are :

1. DPC.py
2. iFeature.py

The github link to the python toolkit: iFeature

# Answer
A few additional hints that may help clarify the algorithm:

The `20` magic number is because there are twenty amino acids. The `400` magic number is because there are `20 * 20` potential dipeptides.

These peptides are then hashed by converting their letter to a number in the range \[0, 20). You can then multiple one of them by 20 and add the other to get a number between \[0, 400). In the case of this function the first amino acid "code" is multiplied by 20. This number can be converted back to the dipeptide sequence when outputting the data. The translation table they use is generated here `diPeptides = [aa1 + aa2 for aa1 in AA for aa2 in AA]`.

The nested loops are to loop through each entry and then each sequence to get every dipeptide sequence, which can be done in linear time.

Just as a quick note on hashing DNA and peptide sequences. You can generally convert them to an array index by using increasing powers of the alphabet size multiplied by the letter index of the alphabet. For example, for DNA you can convert `A=0, C=1, G=2, T=3`, and you can hash fourmers by using `4^4 * a[i] + 4^3 * a[i + 1] + 4^2 * a[i + 2] + 4^1 * a[i + 3]`. In code you may see this performed with bitshifts for DNA since the alphabet size is a power of 2.

The advantage of this hashing scheme is. thateach sequence is guaranteed to have a unique index. The problem with this method is that as if you wish to use longer k-mers the memory usage will grow exponentially. For example, for 16-mers you need 4 billion combinations, which is the size of a 32 bit unsigned integer.
> 3 votes

---
Tags: python, featurecounts
---


# Question
Title: How to plot the gap distribution of contigs wrt to reference genome
I have a contigs file that I generated using Minia and then I have used Minimap2 to map these contigs again to the reference genome. Now I want to plot the gap distribution i.e. gap distance of the set of contigs with respect to the positions in the reference genome. For example: something like , if there are m contigs obtained from the Minia, then X axis will show 0 to m-1 gap points (in an ideal case where there is no overlap between the contigs and the contigs are aligned back to back without space).

The file from minimap2 looks like this:

```
0   17679   0   17679   +   NZ_BCML01000091.1   39426   18736   36415   17679   17679   60  NM:i:0  ms:i:35358  AS:i:35358  nn:i:0  tp:A:P  cm:i:3274   s1:i:17675  s2:i:49 de:f:0  rl:i:75 cg:Z:17679M
1   16011   0   16011   +   NZ_BCML01000021.1   105137  86      16097   16011   16011   60  NM:i:0  ms:i:32022  AS:i:32022  nn:i:0  tp:A:P  cm:i:3008   s1:i:15996  s2:i:0  de:f:0  rl:i:15 cg:Z:16011M
2   2779    0   2779    -   NZ_BCML01000048.1   72397   14902   17681   2779    2779    60  NM:i:0  ms:i:5558   AS:i:5558   nn:i:0  tp:A:P  cm:i:519    s1:i:2751   s2:i:0  de:f:0  rl:i:45 cg:Z:2779M
4   177     0   177     +   NZ_BCML01000150.1   16425   15710   15887   177     177     60  NM:i:0  ms:i:354    AS:i:354    nn:i:0  tp:A:P  cm:i:30     s1:i:160    s2:i:0  de:f:0  rl:i:15 cg:Z:177M

```

where the 1st column is the contig ID and 8th and 9th column show the start and end position of the contig on the reference genome.

I have used the samtools for obtaining a length distribution plot, but not understanding how to do a gap distribution plot. I can use awk on linux or python or R. Anything that is going to work. Thank you for all the help.

# Answer
> 2 votes
## a better solution

I wasn't too happy with the previous answer that I gave, so I was fiddling around, and I realized that `bedtools genomecov` will give you this for free with the following workflow:

```
minimap2 -ax asm10 ref.fasta contigs.fasta | \
samtools sort --output_fmt BAM > contigs_to_ref.bam
# -bga means output in bedgraph format, including all positions
bedtools genomecov -ibam contigs_to_ref.bam -bga > contig_cov.bg

```

If you look at the resulting file, it will report all regions of identical coverage:

```
head contig_cov.bg
# prints the following:
# reference chr/contig, start position, end position, coverage
NODE_1_length_442924_cov_598.199    0   8151    0
NODE_1_length_442924_cov_598.199    8151    8392    1
NODE_1_length_442924_cov_598.199    8392    20203   0
NODE_1_length_442924_cov_598.199    20203   20463   1
NODE_1_length_442924_cov_598.199    20463   44953   0
NODE_1_length_442924_cov_598.199    44953   45216   1
NODE_1_length_442924_cov_598.199    45216   52633   0
NODE_1_length_442924_cov_598.199    52633   52899   1
NODE_1_length_442924_cov_598.199    52899   57698   0
NODE_1_length_442924_cov_598.199    57698   58207   1

```

You can simply filter this file for all regions wherein the fourth column is equal to zero, and take the first 3 columns, and that is a list of gaps in the assembly relative to the reference:

```
awk '$4 == "0" { print $0 }' contig_cov.bg | cut -f1-3
# prints
NODE_1_length_442924_cov_598.199    0   8151
NODE_1_length_442924_cov_598.199    8392    20203
NODE_1_length_442924_cov_598.199    20463   44953
NODE_1_length_442924_cov_598.199    45216   52633
NODE_1_length_442924_cov_598.199    52899   57698
NODE_1_length_442924_cov_598.199    58207   73398
NODE_1_length_442924_cov_598.199    73733   78601
NODE_1_length_442924_cov_598.199    78959   80147
NODE_1_length_442924_cov_598.199    80419   86967
NODE_1_length_442924_cov_598.199    87223   93872

```

You could even do this as a one-liner probably though that's getting a little convoluted.

## Original answer

It's a bit circuitous, but one way to do this would be to follow this answer's suggestion and generate a pileup for the alignment of your contigs to the reference.

You can get minimap2 to give SAM/BAM output with the `-a` flag. You then use the SAM/BAM to generate a pileup or a coverage file, using `samtools mpileup` or `samtools depth -a`. You can then use a heuristic to find contiguous regions with zero coverage of contigs. You could for example loop through a depth file, here is some python code that gives the basic idea:

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

# minimap2 -ax asm10 ref.fasta contigs.fasta | \
# samtools sort --output_fmt BAM > contigs_to_ref.bam
# samtools depth contigs_to_ref.bam > depth.txt
depth_file = "depth.txt"

with open(depth_file) as dep_handle:
  gap_sizes = list()
  start = None
  pos = None
  for line in dep_handle:
    # needs also logic to handle multiple chromosomes in reference
    # but I'm lazy
    fields = line.split()
    if int(fields[2]) == 0:
      pos = int(fields[2])
      if start is None:
        start = int(fields[1])

    elif fields[2] >0:
      if start is None:
        continue
      else:
        gap_sizes.append(pos - start)
        start = None

print(gap_sizes)

```

Probably there is a simpler way to do this but I don't know it.

---
Tags: r, python, fasta, ggplot2, minimap2
---


title: Reproducible PyTorch Model

text:
# Question
Title: Reproducible PyTorch Model
I'm using PyTorch (1.7.1), PyTorch Geometric (1.6.3), NVIDIA Cuda (11.2).

I need to make a neural network reproducible for a competition.

However, when I try:

```
device = torch.device('cuda:0')
rand = 123
torch.manual_seed(rand)
torch.cuda.manual_seed(rand)
torch.cuda.manual_seed_all(rand)
torch.backends.cudnn.deterministic = True
torch.backends.cudnn.benchmark = False
np.random.seed(rand)
random.seed(rand)

```

the results are different each time I run the code. How can I fix it?

(For a reference, it always comes out the same in `device = torch.device('cpu')`.)

# Answer
There is this nice function from Kaggle:

```
def seed_everything(seed=1234):
    random.seed(seed)
    os.environ['PYTHONHASHSEED'] = str(seed)
    np.random.seed(seed)
    torch.manual_seed(seed)
    torch.cuda.manual_seed(seed)
    torch.backends.cudnn.deterministic = True

seed_everything()

```

There are several version of this function for torch, tensorflow and combinations. This specific one is from here.

Run this before you run your model. But not only the fitting function. It needs to be before the validation split function as well.
> 2 votes

---
Tags: machine-learning, deep-learning, pytorch
---


# Question
Title: Site Mutations in a 3D protein structure to locate the core of mutations
I have a pdb file for a protein (3D structure predicted by a tool) and around hundreds of known point mutations in the protein. I want to visualize the protein with all the mutations to see the hub of or core where the mutations are mostly. Is there are tool that does this? I would give the structure in a pdb file and a list of point mutations myself.

# Answer
Clustering mutations could reveal lots about the protein providing one has the functional impact of each mutation (this cluster are benign, this cluster causes misfolding etc). Many tools are emerging to explore this and there are plenty of established tools that can handle this.

VarMap is a tool that is able to map genomic variant coordinates to the best available structure.

PyMol can open PDB files and you can programmatically or manually mutate each position using the mutagenesis tool.
> 0 votes

---
Tags: protein-structure, 3d-structure, point-mutation
---


# Question
Title: Assembling all transcripts for an individual gene? (using single sequence to seed the assembly)
Let's say I have a candidate gene and I believe that in an individual sample, the genome sequence differs from the reference which then interferes with alignment.

Is there a way for me to do a "transcriptome assembly" with a specific read or sequence as a starting point? And only do the k-mer assembly from that single sequence?

I could probably code something myself but I feel like a tool to do this must already exist?

The goal here is to avoid the significant overhead of assembling a whole transcriptome when I'm only interested in one gene.

NOTE: In this case the data is RNA-seq

# Answer
> 1 votes
In essence what you're looking for is targeted assembly. There aren't a lot of tools out there for that that I'm aware of, but on of them is Mapsembler. The general idea is to input a seed region and then let the program perform iterative extension of it. That sort of strategy should largely match your needs.

---
Tags: rna-seq, assembly, k-mer, stringtie, trinity
---


title: mpiblast alternatives?

text:
# Question
Title: mpiblast alternatives?
I am looking for a sequence alignment software capable of pairwise alignment of a large number of protein sequences (10^7).

I tried to go for mpiblast: https://wiki.canterbury.ac.nz/display/UCHPC/Installing+mpiBLAST

but seems their website was down.

Do you know any alternatives to mpiblast which is also highly efficient in protein sequence alignment?

Thanks!

# Answer
The Diamond protein aligner is a very popular tool for pairwise protein alignment which is much faster than blastp while having comparable sensitivity. Diamond also supports the use of using multiple compute nodes (similar to mpiBLAST). You can view the details of the new Diamond 2.0 release in this recent Nature Methods paper (source of the image as well).

> ...we aligned all 281 million protein sequences from the NCBI nr database against the UniRef50 database, which consists of 39 million sequences, using DIAMOND (v2.0.7) in ultra-sensitive mode...
> 4 votes

---
Tags: sequence-alignment, blast
---


# Question
Title: Negative binomial modeling of RNA-Seq data
A common way to model RNA-seq data is using a negative binomial distribution, where each sample-gene pair is modeled by a different negative binomial distribution with mean $\mu\_{ij}$ where $i$ and $j$ are indices for genes and samples, respectively (see article). My question is, given we only have a single observation ($X\_{ij}$) from each sample-gene pair, how come we can learn the $\mu\_{ij}$ from that single observation? We need multiple observations to estimate a mean, right? Maybe I am missing something critical.

# Answer
If reads were independently sampled from a population with given, fixed fractions of genes, the read counts would follow a multinomial distribution, which can be approximated by the Poisson distribution. However, the Poisson distribution assumes the mean equals the variance, which is usually not satisfied by the data. And for the RNA sequencing data, the variance is usually larger than the mean. So we can assume the data follows a $Y|\lambda\sim Poisson(\lambda), \lambda\sim Gamma(r\_0,b\_0)$ hierarchy.

$$\begin{align} P(Y=y)&=\int\_{0}^{\infty}P(Y=y|\lambda)f(\lambda)d\lambda\\ &=\int\_{0}^{\infty}\frac{e^{-\lambda}\lambda^y}{y!}\frac{b\_0^{r\_0}}{\Gamma(r\_0)}\lambda^{r\_0-1}e^{-b\_0\lambda}d\lambda\\ &=\frac{b\_0^{r\_0}}{\Gamma(r\_0)y!}\int\_{0}^{\infty}e^{-\lambda}\lambda^y\lambda^{r\_0-1}e^{-b\_0\lambda}d\lambda\\ &=\frac{b\_0^{r\_0}}{\Gamma(r\_0)y!}\int\_{0}^{\infty}e^{-\lambda(b\_0+1)}\lambda^{y+r\_0-1}d\lambda\\ &=\frac{\Gamma(r\_0+y)}{\Gamma(r\_0)y!}\frac{b\_0^{r\_0}}{(b\_0+1)^{r\_0+y}}\int\_{0}^{\infty}\frac{1}{\Gamma(r\_0+y)}e^{-\lambda(b\_0+1)}\lambda^{y+r\_0-1}(b\_0+1)^{r\_0+y}d\lambda\\ &=\frac{\Gamma(r\_0+y)}{\Gamma(r\_0)y!}\frac{b\_0^{r\_0}}{(b\_0+1)^{r\_0+y}}\int\_{0}^{\infty}\frac{1}{\Gamma(r\_0+y)}e^{-\lambda(b\_0+1)}\[\lambda(b\_0+1)\]^{r\_0+y}\frac{1}{\lambda}d\lambda\\ &=\frac{\Gamma(r\_0+y)}{\Gamma(r\_0)y!}\frac{b\_0^{r\_0}}{(b\_0+1)^{r\_0+y}}\\ &=\frac{\Gamma(r\_0+y)}{\Gamma(r\_0)y!}\left(\frac{1}{b\_0+1}\right)^y\left(\frac{b\_0}{b\_0+1}\right)^{r\_0}\\ \end{align}$$ which is negative binomial distribution.

Negative binomial distribution depicts the number of failures $k$ before $r$th success of a series of independent Bernoulli experiments, with success probability $p$: $$f(k;r,p)\equiv Pr(X=k)=\frac{\Gamma(k+r)}{k!\Gamma(r)}p^r(1-p)^k\quad k=0,1,2,\cdots$$

$$\mu=\mathbb E(X)=\frac{r(1-p)}{p}$$ $$\sigma^2=\mathbb V(X)=\frac{r(1-p)}{p^2}$$

Then $$p=\frac{\mu}{\sigma^2}$$ $$r=\frac{\mu^2}{\sigma^2-\mu}$$ The probability mass function (PMF) is:

$$f(k;r,p)\equiv Pr(X=k)={k+\frac{\mu^2}{\sigma^2-\mu}-1\choose{\frac{\mu^2}{\sigma^2-\mu}-1}}\left(\frac{\sigma^2-\mu}{\sigma^2}\right)^k\left(\frac{\mu}{\sigma^2}\right)^{\frac{\mu^2}{\sigma^2-\mu}}$$

The $\mu$ and $\sigma^2$ of each gene can be estimated from data.

Assume that the number of reads in sample $j$ that are assigned to gene $i$ can be modeled by a negative binomial (NB) distribution,

$$K\_{ij}\sim NB(\mu\_{ij},\sigma^2\_{ij})\tag1$$ which has two parameters, the mean $\mu\_{ij}$ and the variance $\sigma^2\_{ij}$. The read counts $K\_{ij}$ are non-negative integers. $q\_{i, \rho(j)}$ is proportional to the expectation value of the true (but unknown) concentration of fragments from gene $i$ under condition $\rho(j)$. The size factor $s\_j$ represents the coverage, or sampling depth, of library $j$.

* First, the mean parameter $\mu\_{ij}$ , that is, the expectation value of the observed counts for gene $i$ in sample $j$, is the product of a condition-dependent per-gene value $q\_{i, \rho(j)}$(where $\rho(j)$ is the experimental condition of sample $j$) and a size factor $s\_j$,

$$\mu\_{ij} = q\_{i, \rho(j)}s\_j\tag2$$

* Second, the variance $\sigma^2\_{ij}$ is the sum of a shot noise term and a raw variance term,

$$\sigma^2\_{ij}=\mu\_{ij}+s\_j^2v\_{i,\rho(j)}\tag3$$

* Third, we assume that the per-gene raw variance parameter $v\_{i,\rho(j)}$ is a smooth function of $q\_{i, \rho(j)}$,

$$v\_{i,ρ(j)}=v\_ρ(q\_{i, \rho(j)}) \tag4$$

This assumption is needed because the number of replicates is typically too low to get a precise estimate of the variance for gene $i$ from just the data available for this gene. This assumption allows us to pool the data from genes with similar expression strength for the purpose of variance estimation (see article).
> 2 votes

# Answer
The $\mu\_{ij}$ in that section of the DESeq manuscript is the expected value of sample $j$ in gene $i$ given its group association (with expected value $\mu\_{ij}$). This is why $K\_{ij}$ (this is presumably what you meant by $X\_{ij}$) is a function of $\mu\_{ij}$, rather than $\mu\_{ij}$ itself. I think you're just swapping $K\_{ij}$ and $\mu\_{ij}$ in your head.
> 1 votes

---
Tags: rna-seq, deseq2
---


title: creating a dictonary in Python

text:
# Question
Title: creating a dictonary in Python
From the FASTA header I get the following information:

```
gi|1543850138|ref|NC_040192.1| Desmodium mottle virus isolate UG5 segment DNA-B, complete sequence
['gi', '1543850138', 'ref', 'NC_040192.1', ' Desmodium mottle virus isolate UG5 segment DNA-B, complete sequence']

```

Is there a better way to create a dictonary as I did with the below script?

from Bio import SeqIO

organismInfo = {}
for seq_record in SeqIO.parse("data/landplantsviruses.fasta","fasta"):
    print(seq_record.description)
    temp = seq_record.description.split('|')
    print(temp)
    organismID = temp[3]
    organismName = temp[4].lstrip()
    print(organismID, organismName)
    organismInfo[organismID] = organismName
print(organismInfo)
```
</code></pre>

# Answer
You may notice that Uniprot and NCBI have different standard: so your solution is not universal, which is why BioPython does not split the description for you —it could know different 'dialects' of fasta files but that is a tedious minefield.

About coding, one write something in a code golf format, but this is problematic as it is hard to interpret and adapt, i.e. it's a bad practice. Having said that if made a list of keys you could use `zip` with your list of elements as the list of keys may be more intelligible.

```
header = 'gi|1543850138|ref|NC_040192.1| Desmodium mottle virus isolate UG5 segment DNA-B, complete sequence'
keys = 'db|id|type|id2|wordy.'.split('|')
organism_info = dict( zip(keys, map(str.strip, header.split('|')) ))

```

Also, for variables in Python, camel case is for Classes, not instances. Even if you have to send data to JS, using a constistent style, i.e. the PIP guideline one, avoids errors.
> 2 votes

---
Tags: python
---


title: Edit end column in GFF3 file

text:
# Question
Title: Edit end column in GFF3 file
I am trying to edit the "end" of my sequences in a GFF3 that contains 10000 sequences (so I do not want to do it manually). I want to add 30 residues to every single fragment. How can I do this?

Currently, I have a file that looks like this:

```
##gff-version 3
sp_P45452_MMP13_HUMAN   SignalP-5.0 signal_peptide  1   19  0.994300    .   .   .
sp_Q8R2H5_PHLD_RAT      SignalP-5.0 signal_peptide  1   23  0.986089    .   .   .

```

I would like to edit it so that it looks like this:

```
##gff-version 3
sp_P45452_MMP13_HUMAN   SignalP-5.0 signal_peptide  1   49  0.994300    .   .   .
sp_Q8R2H5_PHLD_RAT      SignalP-5.0 signal_peptide  1   53  0.986089    .   .   .

```

Thank you!

# Answer
That's what awk excels at! Simply set the input (`-F '\t'`) and output (`-v OFS='\t'`) record field separators to tab, and add 30 to the value of the fifth field on lines that don't start with a `#`:

```
$ awk -F'\t' -v OFS='\t' '/^[^#]/{$5+=30}1' file.gtf
##gff-version 3
sp_P45452_MMP13_HUMAN   SignalP-5.0 signal_peptide  1   19  30.9943 .   .   .
sp_Q8R2H5_PHLD_RAT  SignalP-5.0 signal_peptide  1   23  30.9861 .   .   .

```
> 3 votes

---
Tags: bedtools, gff3
---


# Question
Title: GTF file format causing GeneMark-ES error
I am predicting genes using `GeneMark-ES`, I managed to generate a `gtf` file, but when I tried to extract the nt/protein sequence using `get_sequence_from_GTF.pl`, I am running into an error in the format of the gtf file, which I find puzzing as it was the gtf file was generated using `gmes_petap.pl` (an script in the `GeneMark-ES` package).

Here is the command I used and the error message which also prints the first line of the genemark.gtf file:

```
$ get_sequence_from_GTF.pl genemark.gtf DT_CAJQTU010000037_1.fasta

error in file format GTF: CAJQTU010000037.1 Desoria tigrina genome assembly, contig: utg000040l, whole genome shotgun sequence     GeneMark.hmm3   gene    2099    4039    .       +       .       gene_id "1_g"; transcript_id "1_t";

```

I think I have found the line in `get_sequence_from_GTF.pl` that is causing the error, but not being a perl monk, I am having troubles to read figure out what is wrong with my `gtf` file and more importantly figure out how to fix it.

```
...
if ( /^(\S+)\t\S+\t(\S+)\t(\d+)\t(\d+)\t.\t([\+-])\t([012.])\tgene_id \"(\S+)\"/ )
    {
    ... (some stuff)
    }
else
{
  if ( /^\s+$/ ) { next };
  if ( /^#/ ) { next };
  die( "error in file format GTF: $_\n" );
}
...

```

# Answer
> 2 votes
I'm afraid this is a clear bug in the tool you are using. The regular expression requires that all fields in the GTF file contain no whitespace. The `\S+` means "one or more non-whitespace characters" but your "source" field, the second field in the GTF file, contains spaces: `Desoria tigrina genome assembly, contig: utg000040l, whole genome shotgun sequence`.

The thing is, that is a perfectly valid field for a GTF file! GTF fields can contain spaces, they just can't contain tabs. I tried to find a repository for the software to suggest a fix, but found only an empty one. In any case, try changing the regex you show to this:

```
#if ( /^(\S+)\t\S+\t(\S+)\t(\d+)\t(\d+)\t.\t([\+-])\t([012.])\tgene_id \"(\S+)\"/ )
if ( /^([^\t]+)\t[^\t]+\t([^\t]+)\t(\d+)\t(\d+)\t.\t([\+-])\t([012.])\tgene_id \"(\S+)\"/ )
    {
    ... (some stuff)
    }

```

I can't guarantee that this will work as expected, please test! The details will depend on what they do with the captured groups (the `some stuff`), but at least that regex will match and I *think* I have kept the numbering of the capture groups.

# Answer
> 1 votes
You might use `agat_sp_extract_sequences.pl` from AGAT instead.

To extract and translate the coding regions:
`agat_sp_extract_sequences.pl -g infile.gxf -f infile.fasta -p`

see here for all possibilities with that script: https://agat.readthedocs.io/en/latest/tools/agat\_sp\_extract\_sequences.html

Otherwise with the modification suggested by @terdon, your script should work fine.

---
Tags: gtf
---


# Question
Title: How do I call variants using BCFTools without piping mpileup file?
I need to call variants on a large number of reads using both `VarScan2` and `BCFTools`. Both of these variant calling tools use `mpileup` file as input to call the variants. `VarScan2` command is straightforward, as we can specify the `mpileup` file as a command parameter. The problem is, `BCFTools` expects `mpileup` to be piped from another command and I cannot find an input parameter to specify a `mpileup` file that already exists (the one used with `VarScan2`). It seems horribly inefficient to run the same command twice instead of reusing the file obtained from previous step (especially because generating `mpileup` takes a substantial amount of time on my computer).

Is there a way to reuse an existing `mpileup` file with `BCFTools`?

I am using `Nextflow` to run my pre-processing pipeline, so could it be possible that I can pipe an input channel with `mpileup` into the `BCFTools`'s `call` command?

~~I also tried doing `cat file.mpileup | bcftools call ...` but this does not seem to work.~~

Edit: `cat file.mpileup | bcftools call ...` does work, which would be a solution to this question, however, answers offer other, even more elegant, solutions.

# Answer
Why not use `tee`? Please see this post.

You can write `mpileup` to file while also piping it to another process. For example, `samtools mpileup ... | tee file.mpileup | bcftools call ...`

the `bcftools` will run while generating the `file.mpileup`. After it's done, you can use `file.mpileup` for your `VarScan2`.

You can test it with the following:

```
touch helloworld.txt
ls | tee file.txt | grep helloworld    # output helloworld.txt
cat file.txt | grep helloworld         # output helloworld.txt

```
> 2 votes

# Answer
Bcftools and samtools both follow the very common Unix filter pattern: they read from all input files specified (as non-option plain arguments), or if there are none they read from standard input.

`bcftools call` is no exception:

```
$ bcftools call

About:   SNP/indel variant calling from VCF/BCF.
         To be used in conjunction with bcftools mpileup.
         [etc]
Usage:   bcftools call [options] <in.vcf.gz>

```

So you can write either of the following equivalently:

```
cat file.mpileup | bcftools call [options]
bcftools call [options] file.mpileup

```

This assumes that the file is VCF/BCF as produced by `bcftools mpileup`, and not actually textual “mpileup output” as suggested by the filename. If it is the latter that VarScan2 is expecting, you won't be able to share the input files between the two after all.

See this workflow example, which demonstrates `bcftools call` both with a pipe and a previously made input file.
> 1 votes

---
Tags: variant-calling, nextflow, bcftools, mpileup
---


# Question
Title: Is it possible to do DEG analysis without replicates?
I have three bulk RNA-seq samples (test1, test2, test3) without replicates.

And I noticed that DEG analysis tools such as DESeq2/edgeR/etc cannot be applied for data with no replicates.

So, I just drew a heatmap with expression value (TPM) of genes of interest, using NMF (aheatmap).

But now I want to perform DEG analysis across the three samples. Is there anything I can try for the analysis?

My TPM data (.csv) looks like below.

Thank you so much.

# Answer
According to DESeq2, you cannot analyze your data without having replicates.

In short, you have an experiment where you compare A to B because, in your theory, A is more expressed than B. However, you do the experiment only once. Without repeating the same experiment again under the same circumstances, there's no way to confirm you will have the same results and that your theory is correct.

You can still compare the ratio of `A/B` but that will not be meaningful in any way since there's not enough statistical power behind it. In your case, doing the ratio of `test1/test2`, `test2/test3` and `test1/test3`.
> 2 votes

# Answer
In this situation, you can apply the strategy in section 2.12 in the edgeR manual which basically "makes up" a dispersion estimate and then run the normal DE strategy. This is obviously neither reliable, nor publishable but at least it gives you a list of genes you can use for validation. Treat results with care, focus on genes with decent counts (not low on both conditions of the contrast) and select those with large effect sizes (fod changes). Alternatively, you can run `vst` from DESeq2 and then rank them by fold change, again focusing on genes with decent (non-low) counts. Both is not really reliable, so be careful with the results. As you have no information on the variability of the groups the fold changes are of limited value, but if you are forced to go with these data it is better than nothing.
> 3 votes

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


# Question
Title: Defining a bedfile for chromosome X without pseudo-autosomal regions
I would like to define a bedfile for chromosome X, to target only basepairs outside of pseudo-autosomal regions (and going forward other regions, too).

First, I wasn't able to find an already prepared bedfile in any of the online repositories (but let me know if it is available somewhere).

Second, I compiled the bedfile for all basepairs in chromosome X as a starting point (by looking at the latest assembly at genome browser):

```
chrX    0       156040895       # as defined at https://genome.ucsc.edu/

```

The genome browser is 1-based and the last base is inclusive, while the bedfile is 0-based and last base is not inclusive, therefore I had to subtract 1 from the bed `start`.

Afterwards, I wanted to exclude the pseudo-autosomal regions as defined at this page. However, I am not sure if basepair "Start" and "Stop" are 0- or 1-based and inclusive/not inclusive. I might use the same convention used on the genome browser, but I am not sure if that's correct in this case. I wasn't able to find this information on their website or online.

Assuming that the NCBI "Start" and "Stop" field follow the UCSC parameters (1-based and inclusive), the resulting bed (PAR regions excluded) might be the following:

```
chrX    0       10000
chrX    2781479 155701382

```

**Edit**: the section on the centromeres has the "Length" parameter, from which I could infer that "Start" and "Stop" base-pairs are both inclusive. At the same time, the "Stop" base-pair of PAR#2 in chromosome X is base-pair 156,030,895, which is also the last base-pair number if we were considering the UCSC genome browser notation (1-based and inclusive). Therefore, it seems to me clear that the same notation is used, and the proposed bed-file should be correct.

# Answer
> 2 votes
The UCSC Genome Browser renders data visually with a one-based, closed index. The underlying data are often stored and presented as files worked with on the command-line, using zero-based, half-open coordinates.

This difference in indexing can be important when using the website to view data or when using the download tools to grab a dataset for offline calculations, such as files in BED format.

There is a good overview of indexing on Biostars.

To solve this in a general way:

The UCSC Kent Utilities `fetchChromSizes` tool can be used to retrieve the sizes of chromosomes for your genome assembly of interest, e.g. `hg38`.

Using this assembly, to build a sorted BED file for set operations:

```
$ fetchChromSizes hg38 | grep -v '_*_' | awk -v FS="\t" -v OFS="\t" '{ print $1, "0", $2 }' | sort-bed - > hg38.bed

```

Replace `hg38` with your assembly of interest.

To then filter this for `chrX`:

```
$ awk '($1=="chrX")' hg38.bed > hg38.chrX.bed

```

The intervals in `hg38.bed` and `hg38.chrX.bed` use a zero-based, half-open coordinate system.

To then get regions of `chrX` without autosomal regions, you could define a sorted BED file containing those regions, which also use zero-based, half-open coordinates:

```
$ echo -e 'chrX\t10000\t2781479\nchrX\t155701382\t156030895' > pseudoautosomal.chrX.bed

```

I am assuming that the coordinates on Ensembl and NCBI use one-based indexing. Converting to a zero-based index requires subtracting one base from the start positions of these intervals.

It would be ideal if these groups were to clarify the indexing they use. I am making this assumption based on how Ensembl defines regions on `chrY`, which start at position `1` and not `0`, as one would expect to see with a one-based coordinate system.

The numbers reported on Ensembl and NCBI are otherwise identical, so I also assume that NCBI's numbers use the same indexing and require the same correction to the start base.

In any case, to calculate the difference between these regions and `chrX`:

```
$ bedops --difference hg38.chrX.bed pseudoautosomal.chrX.bed > answer.bed

```

This is what my output looks like:

```
$ more answer.bed
chrX    0   10000
chrX    2781479 155701382
chrX    156030895   156040895

```

These coordinates are in BED format and so use a zero-based, half-open coordinate system, by convention.

To turn these into one-based, closed coordinates — to query signals over exact regions in the UCSC Genome Browser, for instance — you would add one base to the start positions of each genomic interval.

(Edit: Fixed PAR#2 start coordinate per comment.)

---
Tags: bed, chromosomes
---


# Question
Title: Why does the SARS-Cov2 genome has letter t
`ATTAAAGGTT TATACCTTCC CAGGTAACAA ACCAACCAAC TTTCGAT...` is part of the 5'UTR of genome of an RNA virus SARS-Cov-2.

RNA contains letters `C`+`G`+`A`+`U`.
DNA contains letters `C`+`G`+`A`+`T`.

Why are there letters `T`, in the SARS-Cov-2 genome, instead of `U`?

# Answer
We sequence and therefore typically report assemblies as DNA sequences, even if they're actually RNA.
> 10 votes

---
Tags: genome, phylogenetics, rna, sars-cov-2, covid-19
---


# Question
Title: How to calculate coverage depth of genes in an annotated sequence?
I have a metagenomic dataset across multiple samples that has been grouped and binned (we did this due to poor sequence quallity).

It is several samples taken from termite guts of differing castes, the samples have been grouped to make binning easier due to poor sequence suitability. They have been binned via ASVs, assembling reads from all samples to generate the putative bacterial species. So I have fasta files of putative species that have been annotated via prodigal. I want to now calculate the abundance of each bin in a given sample. I think the best way to do this would be at the gene level by looking at the number of reads that map to a gene in any given sample. I am just unsure how to do this.

I have annotated the bins and now want to calculate the abundance of each bin in a sample but looking at the number of reads from a sample that map to each gene in a bin.

Is this a viable method? If so how would I do it and if not what else could I do?

# Answer
The best way to approach it would be to take trimmed reads map them to the bins (i.e. with bowtie2) and estimate the abundance of each bin with samtools by getting the depth of coverage:

https://www.metagenomics.wiki/tools/samtools/breadth-of-coverage

or you can use QUANT\_BIN module from the METAWRAP to do it for you. Follow this for the METAWRAP approach:

https://github.com/bxlab/metaWRAP/issues/40
> 1 votes

---
Tags: genome, metagenome
---


title: Generating sequences from HMM

text:
# Question
Title: Generating sequences from HMM
I encountered this question in a textbook (source: https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.721.2540&rep=rep1&type=pdf) and was trying to workout the solution for it. Was wondering for the last question, wouldn't the probability be 0?

My working:

START -\> M1 -\> M2 -\> M3 -\> M4 -\> M5 -\> M6 -\> M7 -\> END

Probability Vector \[A, T, C, G\]

M1 : \[0.6, 0.4, 0, 0\]

M2 : \[0.6, 0.4, 0, 0\]

M3 : \[0.4, 0.4, 0.2, 0\]

M4 : \[0.6, 0.4, 0, 0\]

M5 : \[0.4, 0.6, 0, 0\]

M6 : \[1, 0, 0, 0\]

M7 : \[0, 0, 0, 1\]

To generate CTAAAAG = 0 * 0.4 * 0.4 * 0.6 * 0.4 * 1 * 1 = 0

# Answer
> 1 votes
*\[Comments from other post migrated as answer\]*

The parenthetical statement in the first bullet point says **without corrections**. In which case, if an unseen nucleotide has a probability of zero, the probability of the sequence is zero. If the unseen nucleotides had the probability of NaN, then the sequence has a probability of NaN, which is a more thorough solution.

The correction would be adding a small iota value to the probabilities to be non zero —**fudging** against the zero outcome you otherwise find, which may be the result of the training being too small. You have 5 training sequences, so the iota would be \<1/5. Let's hypothetically assume that if you were to have picked a 6th MSA sequence it would have had the missing value, so you could set the iota at 1/6. The assumption that you may find that element were you to sample twice as many is more statistically sensible, so iota to 1/10.

I used the letter iota, which is common parlance for something small, but iota is a really bad in maths as it looks like a 1 with English handwriting or in print: you can call the fudge value anything.

---
Tags: homework, hidden-markov-models, multiple-sequence-alignment
---


title: Variant vs Allele vs SNP

text:
# Question
Title: Variant vs Allele vs SNP
Coming from a CS background. Reading through the wikipedia page, these all sounds like the **same** thing: Variant, Allele, and SNP.

> **Variant/Allele/SNP:** Some gene locus that differs from the idea human.

For example, if 99% of humans have a T at some locus, and you have a G, C, or A, then you have a Variant/Allele/SNP at that locus.

Is there intuitive way to explain the difference?

# Answer
Variant = a locus in the genome where there are differences between individuals

Allele = one of the possible bases/sequences that can occur at the variant

SNP = single nucleotide polymorphism, a type of variant. Here the alleles are just single bases, but other variants can be insertions/deletions of one or more bases.
> 4 votes

# Answer
First some basic biology, taking human as an example. Humans have 22 chromosomes (image source):

As you can see (if you squint) in the image, we have two copies of each of the 22. With the exception of the sex chromosomes, X and Y (seen in the bottom right of the image above), you have the same complement of genes on each of the copies of the chromosome. Traditionally, the term "allele" refers to the sequence of one of the two copies.

We inherit one copy of each chromosome from each of our parents. So, if you have inherited a mutation at position $Ψ$ of chromosome 21 from one of your parents, you will be "heterozygous" for that mutation. This means that of the two copies of chromosome 21 in your cells, one will have the mutant *allele* at position $Ψ$, and the other the wild-type (i.e. "normal", although that's possibly a misnomer) *allele*.

SNP stands for "single nucleotide polymorphism". This means that it describes a difference *of a single nucleotide* between the genome you are looking at and the reference genome. For example, if the reference genome has a G at position Ψ of chromosome 21 and you have a T, then this is a SNP.

The term SNP is also commonly used to refer to polymorphic positions. This means specific, single base positions of the genome where variation is common in the human population. Having differences from the reference genome in such polymorphic sites is usually not problematic. The specific use of the term varies a little between subfields. If you're working on population genomics, you are likely to use SNP only for such known polymorphic sites, while others may use *SNP* to refer to any single nucleotide variant.

Which brings us to *variant*. This is a blanket term encompassing any and all forms of variation in a genome. SNPs are variants, but so are larger insertions or deletions (indels), copy number variants, large transposed sequences, the list goes on. You can basically think of *variant* as the catch all term used when you don't want or need to use anything more specific.

To sum up, we have:

* Allele: one of the two copies of a genomic sequence present in (diploid) cells.
* Variant: a site where the genome you are looking at differs from a reference genome of that species.
* SNP: a single nucleotide variant often, but not always, in a position of the genome where a significant proportion of the species' individuals show variation.
> 2 votes

# Answer
Just to add on to previous answers, it is important that allele is one of the possible versions or sequences that can happen at a given position.

Take SNP. Most SNPs have two alleles, meaning that in the population we only see two possible nucleotides at this position. For instance, everybody either has an A or a G at this position. This type of SNPs are called **biallelic**.

In rare cases, a SNP can have three different alleles, like A/G/C. These are called **triallelic**. So already you can see that **allele** is different from **SNP**. They are actually completely different concepts.

Also, the word allele has been used for a long time. Scientists first used this to refer to a certain version of genetic variation at one place in a genome. This happened even before we know what is DNA. On the contrary, SNP was only discovered after we can sequence human DNA, which is basically in the past 20-30 years. Therefore sometimes it is hard to distinguish the subtle differences between these terms, because they can carry both original meaning based on primitive understand of heredity as well as newer meaning in the era molecular genetics. It might be helpful to think about what these term are referring to in different time and contexts.
> 1 votes

---
Tags: genetics
---


# Question
Title: Why do SAM and BAM use different coordinate systems?
BAM files are, at least as far as I know, simply binary compressed versions of SAM files. They have the exact same information and are used in the same way. Why then does the SAM format use a 1-based coordinate system while BAM uses a 0-based one?

This is the relevant section of the SAM format specification (emphasis mine):

> **1-based coordinate system** A coordinate system where the first base of a sequence is one. In this coordinate system, a region is specified by a closed interval. For example, the region between the 3rd and the 7th bases inclusive is \[3, 7\]. **The SAM, VCF, GFF and Wiggle formats are using the 1-based coordinate system**.

> **0-based coordinate system** A coordinate system where the first base of a sequence is zero. In this coordinate system, a region is specified by a half-closed-half-open interval. For example, the region between the 3rd and the 7th bases inclusive is \[2, 7). **The BAM, BCFv2, BED, and PSL formats are using the 0-based coordinate system**.

The two coordinate systems are a very common source of error in bioinformatics, so there must be some reason why these two largely equivalent formats use different systems. For that matter, why do the VCF and BCF formats, whose relationship is equivalent to that between SAM and BAM, also do this? I am guessing there is a good reason and it is equally applicable to both pairs of sister formats.

Note that the `samtools` manual states that "all position coordinates are 1-based" when giving regions, so it must handle the conversion internally. Indeed, I can confirm that this command returns identical md5sums for both inputs:

```
md5sum <(samtools view -b sample.bam chr1:36931633-36931871 | samtools view) \
       <(samtools view sample.bam chr1:36931633-36931871)

```

# Answer
Arithmetic on a zero-based coordinate system is less complicated than that on a one-based system, so it appears zero-based is often (not exclusively) used for binary data formats, like BAM or bigBed, or text formats like BED, where computers are used to more efficiently calculate lengths of or set operations on intervals.

A more complete answer on the benefits of zero-based accounting is available here, which references the original note (EWD831) on this subject from Dijkstra.

When reading other textual data formats like SAM or VCF, a one-based index is perhaps more "natural" for humans used to counting digits on their hands. If you visually inspect an annotation on the closed interval `[3,7]` on a line that denotes some abstract chromosome and starts at position `1`, it is clear to you where those data points are positioned. Colloquially, you understand that the interval starts on the *third* position and ends on the *seventh*. Or, perhaps, people would typically label the initial item in this set the *first* item, and not the *zeroth* item.

A one-based index therefore seems to get used often with formats used for visual inspection or presentation of genomic intervals or sequences, such as in genome browsers — or via CLI tools like `samtools`, when exporting reads to human-readable SAM.

The language used for working with and organizing data can determine choice of indexing, as well. For programs written in C and offshoot languages, including an API like `htslib`, binary data is typically worked with as chunks of memory accessed via pointers or arrays, and it is convention to iterate over that data in increments that use a zero-based index. Contrariwise, R's roots in FORTRAN and focus on linear algebraic operations (and matrix indexing) may perhaps explain why that language uses one-based indices for its lists, matrices, and vectors.

The designer of the BAM and SAM formats, Heng Li, comments below(\*):

> I think 1-based is easier for biologists and mathematicians while 0-based is easier for C/C++/Java programmers. SAM/VCF is for biologists to read. BAM/BCF is for a small number of programmers. I wrote the initial specs of sam, bam and bcf2. PS: when sam and vcf are stored in memory by htslib, everything is 0-based in accordance with C conventions. bam and bcf are largely direct dumps of their in-memory representations.

*Note: Original comments were unfortunately deleted from answer by someone, but the quotes were taken verbatim at the time they were written.*
> 12 votes

---
Tags: sam, file-formats
---


# Question
Title: Is there any database about the abundance of proteins in urine of healthy patients?
I have a dataset about adundance of urine proteins from cancer patients. I would like to compare these abundance levels to healthy patients' urine proteins. However I don't have any data with healthy patients. Is there any database about it?

# Answer
If you have access to the UK Biobank dataset, then they have urine data for healthy individuals.

https://www.nature.com/articles/s41588-020-00757-z
> 3 votes

---
Tags: proteins, database, cancer
---


title: Help with Limma-model

text:
# Question
Title: Help with Limma-model
I'm starting to use Limma to find differential expressed sites within my data and would like to ask if my approach was correct.

The data has a control and a treatment set, each containing data of 48 patients, which were measured before treatment (t0) and after treatment (t1). So I have Control-t0 and Control-t1 as well as Treatment-t0 and Treatment-t1, whereas Control-t0 and Treatment-t0 should not show any effect. Now I would like to know if there is a significant change in Treatment compared to Control. Therefore, I think, I have to adjust for the baseline (t0), right?

The code below is probably the most verbose option Limma provides, but I wanted to use the long method first before I use the shortcuts. So would the Interaction-coefficient I defined within my cont-matrix be the difference between Control and Treatment adjusted for the baseline? Also, I had to turn of adjusting for multiple testing here as my example data is random

```

#Data.frame with columns representing samples and rows different sites/genes
data_df <- data.frame(replicate(8, sample(100,1000, rep=TRUE)))
colnames(data_df) <- c("Patient_01_t0", "Patient_01_t1",
                       "Patient_02_t0", "Patient_02_t1",
                       "Patient_03_t0", "Patient_03_t1",
                       "Patient_04_t0", "Patient_04_t1")

#Metadata for data_df
meta <- data.frame(Sample= colnames(data_df),
                   ID= c(rep("Patient_01",2),
                         rep("Patient_02",2),
                         rep("Patient_03",2),
                         rep("Patient_04",2)),
                   Time= c(rep(c("t0","t1"),4)),
                   Group=c(rep("Control",4),
                           rep("Treatment",4))
)

#Create levels for model-matrix
time <- factor(meta$Time, labels=c('t0','t1') )
group<- factor(meta$Group, labels=c('Control','Treatment') )
lev=paste(group,".",time,sep="")
lev

design <- model.matrix(~0+lev)
design

#Create contrast matrix with all factors of interest
cont_matrix <- makeContrasts( t1_VS_t0_Ctrl = levControl.t1 - levControl.t0,
                              t1_VS_t0_Treat = levTreatment.t1 - levTreatment.t0,
                              Ctrl_vs_Treat_at_t0= levControl.t0 - levTreatment.t0,
                              Ctrl_vs_Treat_at_t1= levControl.t1 - levTreatment.t1,
                              Interaction = (levControl.t1-levControl.t0)-(levTreatment.t1-levTreatment.t0),
                              levels=design)

#Calculate different expressed sites
fit <- lmFit(data_df, design)
fit <- contrasts.fit(fit, cont_matrix)
fit.eBayes <- eBayes(fit, robust=TRUE)

summary(decideTests(fit.eBayes, adjust.method="none"))

```

# Answer
> 3 votes
It looks good to me. One think I could suggest is to see the distribution of your data/model (with limma/voom).

After you get the list of your genes with:

```
top.table <- topTable(fit.eBayes, sort.by = "P", n = Inf)

```

Maybe you can check your multidimensional scaling (MDS):

```
plotMDS(your_data, col = as.numeric(group))

```

Then you can also check with voom

```
y <- voom(data, model.matrix, plot = T)

```

I hope this helps

---
Tags: r, differential-expression, limma
---