# Question
Title: Best data repository to publish a large 'plotly' table containing all annotations on a transcriptome?
I am not sure whether this is the best SE community for this question.

I want to publish a transcriptome paper, along with interactive materials enabling readers to peruse the data behind the discussion. The R package 'plotly' enables rendering online-publishable tables that can be filtered & searched for specific entries.

I would like to publish such a table as supplementary to a paper (e.g. preprint version) however usually data repositories allow for uploading files and text.

Would anyone here know an official data repository with a format such that would enable publishing online an interactive table?

A good example of such a table is given here.

Thanks in advance.

# Answer
You may try Code Ocean. Currently, it is not that easy to use IMHO, but one of their tech leads promises to substantially improve it towards the year end.
> 3 votes

# Answer
I have seen some (1 or more) articles with interactive data displays in the F1000Research journal I don't know recall now if they had interactive tables or plots but they used plotly for sure.
> 2 votes

# Answer
The example looks a lot like a shiny application. Your best bet then is one of the following:

* Ask your IT department to setup a shiny server for you. They can use the free version.
* Post your shiny app on shinyapps.io. Depending on how popular the page becomes, you might have to pay for an account.

The free version of shiny works for most cases (it's what we use), the only problem with it is that it doesn't scale (but if you care about that then you need to either pay for hosting on a pro server or buy a shiny pro license).
> 1 votes

---
Tags: r, public-databases, public-dataset
---


# Question
Title: Article GEO gives counts that are not integers. Should UMI counts be integers?
I thought that `UMI` counts are always integers, but when I opened several datasets provided by `GEO` I got confused because in some articles they provide `counts` that are supposedly non-normalized raw `UMI` count values, but they are not integers. For instance: GSM2043231 or GSE77067

Could anyone clarify, what these people mean by `abundance measurements`?

# Answer
> 3 votes
The data processing section doesn't mention using UMIs at all, so it's likely there were none (whether they'll be present depends on how the libraries are made, with lower throughput methods being less likely to include them). The values are stated as being expected counts from RSEM, which is why they're not always integers.

---
Tags: scrnaseq, umi, geo
---


# Question
Title: Output file names are not known beforehand
I have a folder where the `hdf5` files inputs are (matrices of `scRNAseq` raw `UMI` counts). There is a `python` script that loads them, splits them into multiple matrices, normalizes them, and saves as `csv` files. When the name pattern of the output files is known, things are easy, I know how to do that, but when one input (`hdf5` file) can produce various number of outputs (`csv` files) I am not sure how to approach the issue. Currently my `Snakemake` file relevant part looks the following way:

```
hdf5_files,=glob_wildcards('/main/{hdf5_file}.hdf5')

rule all:
  input:
--->  expand('data_files/normalized/{hdf5_file}.csv', hdf5_file=hdf5_files)
    ...

rule normalization:
  input:
     script = 'python/normalize.py',
     path = '/main/{hdf5_file}.hdf5'
  output:
--->  path = 'data_files/normalized/{hdf5_file}.csv'
  shell:
     "python {input.script} -data {input.path} -path {output.path}"

```

Above I highlighted the `2` lines that I am not sure now how I should write. Any suggestions would be greatly appreciated.

# Answer
Could you give an example of names of input files and output files of `normalize.py` ?

It seems to me that your case could be similar to the one described about the `dynamic` output:

https://snakemake.readthedocs.io/en/stable/snakefiles/rules.html#dynamic-files

Here you would probably use it for the output of `normalization` and for the input of `all` (with a wildcard different from `hdf5_file`).

Disclaimer: I haven't had the opportunity to use this feature so far.
> 4 votes

---
Tags: snakemake
---


# Question
Title: How to map short sequences to long reads, recovering all multiply-mapped high-quality matches
### The dilemma:

I have around one million short sequences (21 bp to several 100s of basepairs) for which I need to identify all occurrences of in 20-30x coverage noisy long reads (both pacbio and ONT). All of the short sequences and long reads are derived from the same individual, and there are enough short sequences given the organism's genome size that each long read should have at least one or multiple short sequences.

### An example:

Given three short sequences, `AA` & `BBBB` & `CCCCCC` , how to identify all occurrences of these sequences in long, uncorrected reads.

```
Sequences to identify in long reads:
  - AA
  - BBBB
  - CCCCCC

Long Reads:
  - 1: ------AA-----------------------CCCCCC------
  - 2: ------AA--------BBBB-----------------------
  - 3: ---BBBB-----------CCCC--
  - 4: ------------------AA-----

```

In this example, the output sam file will contain high MAPQ hits for short sequence `AA` to reads (1, 2, 3), for short sequence `BBBB` to reads (2, 3). For short sequence `CCCCCC` there will be a high MAPQ hit to read 1 and a lower MAPQ hit to read 3

### Limitations:

* Due to the requirements of this project, I cannot perform a de novo assembly then map the short sequences to the assembly output. It is necessary to directly map the short sequences to the long reads.
* I cannot correct the long reads before I map/align the short sequences to them.
* The alignment/mapping should sensitively detect the short sequences in the long reads.
* The short sequences may be anywhere from 1Mbp to 30Mbp in total bases.
* The long reads may be anywhere from 3Gbp to 90Gbp in total bases. So the mapper/alignment technique is hopefully fast. Sensitivity is more important though.

### Research so far:

* `bwa mem` does not output multiply-mapped in a predictable manner that fits this use case https://www.biostars.org/p/304614/
* `blat` seems to output many hits for a single short sequence, but there are no out-of-the-box parameters for mapping short reads to long reads.
* Maybe multiple sequence alignment is the most sensitive, but wouldn't that entail running (no. of short seqs) * (no. of long reads) alignments?

# Answer
> 3 votes
LAST has given the best results for me when I've tried to do this, although I agree with @user172818 that it's not a good idea to map really short reads. This is due to a combination of natural sequence duplication in long DNA sequences (e.g. see here), as well as abundant base calling differences present in single-molecule sequencing.

Minimising error is not necessarily going to be the best option, and concentrating only on unique hits will miss a lot of signal. There are frequently multiple identically plausible positions, even at zero error, and the long reads will have their own associated errors.

I also find the limitations a little odd. Canu can do read-level correction on long reads based on overlapping reads, and if assembly *is* possible, then nanopolish can be used to correct nanopore reads for systematic base-calling error introduced by the methylation of unamplified template.

---
Tags: sequence-alignment, software-recommendation, phylogenetics, long-reads, homework
---


# Question
Title: How to combine multiple kernels of large sample datasets?
I have multiple large sample datasets in matrix format (each has 15000 rows and 5-50 columns) corresponding to different experiments. Each matrix contains the same number of samples(rows) but the variables(columns) are not the same. My objective is to cluster the samples on the basis of all the experiments.

I tried to use Unsupervised multiple kernel learning (UMKL) to integrate the datasets followed by kernelPCA using "mixkernel" package in R (https://www.ncbi.nlm.nih.gov/pubmed/29077792). The UMKL step calculates kernels for each dataset and combines the kernels using 3 different approaches: 1) calculating a consensus kernel from multiple kernels 2) calculating a sparse kernel preserving the original topology of the data 3) calculating a full kernel preserving the original topology of the data

The kernel calculation step was fine but the kernel integration step (all 3 approaches) runs very long and my computer hangs.

Is there any way to handle this problem? More specifically, is there any way to handle multiple kernel integration for large sample datasets?

Any suggestion alternative to using kernel methods will also work.

# Answer
Welcome to Bioinformatics Stackexchange @SD1024. A recently published algorithm known as MOFA (Multi-Omics Factor Analysis, paper, github) is generating a lot of interest, and is designed answer exactly the sort of biological question you are describing. It claims to extract axes from multiple matrices with overlapping samples, but not necessarily overlapping features and handle missing data elegantly. I havn't tried it myself you, but I'm quite excited to do so.
> 1 votes

# Answer
You can use regularized general canonical correlation analysis to cluster those matrices. I would suggest RGCCA, which allows to do this for an arbitrary number of matrices but it doesn't use kernels. There has been some efforts to use kernels on the canonical correlation analysis but there isn't any implementation (yet), you can read more on the related literature, specially the Kernel Generalized Canonical Correlation Analysis.

However in this kind of analysis you will need to come up with a design of how those matrices are related between them. It is not like a PCA where you give some data and you get some data back that summarize the matrix.
> 0 votes

---
Tags: r, clustering, multi-omics
---


title: python package for NNI neighbors

text:
# Question
Title: python package for NNI neighbors
I am working on protein sequence data files for reconstructing phylogenetic tree and I need to generate all **NNI-neighbours** of the tree (two trees are NNI-neighbours if one can be transformed into another by one nearest neighbour interchange operation). Is there any python **package/sorce code** I can access to, which can retrieve all NNI neighbors of a tree. Earlier I have used biopython 1.72 to generate tree from my dataset. But for now I'm looking for a python package which can do NNI search on my tree and return all NNI-neighbors.

# Answer
Both the ETE Toolkit and DendroPy offer methods of node removal and reattachment - I would think it would be easy enough to hack together an NNI movement function with them.

You could also investigate the source code of BioPython's `NNITreeSearcher()` function.
> 2 votes

---
Tags: python, phylogenetics
---


# Question
Title: Getting all NNI trees of a parsimony tree
I have an aligned protein sequence file which I have been using for reconstructing a parsimonious tree. I am currently using NNITreeSearcher.\_get\_neighbors method from Biopython 1.72 but it's way to find one best scored tree only.

```
def _nni(self, starting_tree, alignment):
    """Search for the best parsimony tree using the NNI algorithm (PRIVATE)."""
    best_tree = starting_tree
    while True:
        best_score = self.scorer.get_score(best_tree, alignment)
        temp = best_score
        for t in self._get_neighbors(best_tree):
            score = self.scorer.get_score(t, alignment)
            if score < best_score:
                best_score = score
                best_tree = t
        # stop if no smaller score exist
        if best_score >= temp:
            break
    return best_tree

```

I need to generate all NNI-neighbours of the trees (two trees are NNI-neighbours if one can be transformed into another by one nearest neighbour interchange operation).

1. I want to know how can I get total number of trees generated before the best tree is reached.
2. Is it possible to return all the trees generated before the best tree is reached, so that the all/required number of trees(other than the best tree) may be printed?

# Answer
> 4 votes
As I've studies the `'TreeConstruction'` module of Biopython, I came to know that `'_get_neighbors'` method is the one that is used to get all neighbor trees of the given tree. This is all I've concluded from my findings. If anyone wants to add more, that'll be positively appreciated.

---
Tags: python, phylogenetics, algorithms
---


# Question
Title: filter out FASTQ reads which are shorter
I tried to filter out FASTQ reads which are shorter than 259 bp with

```
bioawk -cfastx 'length() >= 259 {print "@" " " "\n""\n+\n"}' good/SZ005_NoIndex_L002_R1_009.good.fq.gz \
  | gzip >  good-filtered/SZ005_NoIndex_L002_R1_009.good-filtered.fq.gz`.

```

Unfortunately, the output was wrong:

```
@

+

@

+

@

```

What did I miss?

Thank you in advance,

# Answer
> 0 votes
If you have paired-end reads, this solution keeps the two files in-sync (i.e. discard pairs where one of the two reads is shorter than 259). Also, it uses only Unix tools without requiring external programs:

```
paste <(zcat S1_L001_R1_001.fastq.gz | paste - - - -) <(zcat S1_L001_R2_001.fastq.gz | paste - - - -) \
| awk -v FS='\t' 'length($2) >= 259 && length($6) >= 259' \
| tee >(cut -f 1-4 | tr '\t' '\n' | gzip > r1.fq.gz) \
| cut -f 5-8 | tr '\t' '\n' | gzip > r2.fq.gz

```

# Answer
> 3 votes
You need to use variable `$seq`. For example:

```
bioawk -cfastx 'length($seq)>=259{print "@"$name"\n"$seq"\n+\n"$qual}' test.fq.gz

```

There are also more convenient tools like seqtk and seqkit. With those, you may

```
seqtk seq -L260 test.fq.gz

```

---
Tags: fastq, bioawk
---


# Question
Title: Curl error when installing Bioconductor packages with Microsoft R Open 3.5.0
I'm trying to install Bioconductor packages with Microsoft R Open 3.5.0 seems to fail not sure what is the issue.

This is what I'm trying to do

```
source("https://bioconductor.org/biocLite.R")
biocLite("BiocInstaller")

```

But I keep getting an error about "curl" had non-zero exit. How can I install the packages of Bioconductor?

# Answer
I can't comment yet but from what you said it seems that you need to install additional libraries like the one you mentioned: "curl".

from the R terminal you should type:

```
install.packages("curl")

```

then R will install the package. Not sure if there are some other dependencies. If yes R will return you some error message in which there are some package names you need to install.

EDIT:

of course, after installing all the required packages you should be able to install biocLite too.
> 2 votes

---
Tags: r, bioconductor
---


# Question
Title: How to plot circular chromosomes with Gviz?
I am using Gviz to plot several tracks. I create a `GenomeAxisTrack()` and then several `AnnotationTrack()`s that I feed to `plotTracks()`. This produces a stacked output of tracks with annotated arrows as I want but the displays are linear. In some cases I have the complete genomes and I know they are circular (like some bacterial and viral chromosomes), how can I display the plots as circular?

I tried feeding to `AnnotationTrack()`; a `GRanges()` object with a `Seqinfo()` object specifying the `isCircular=TRUE` flag, but this didn't change the resulting plot.

Here is minimal code to create a linear plot, I would like one where the first and last arrow join up in a circle:

```
library(Gviz)
library(stringr)

track <- AnnotationTrack(start=c(1,5,7), end=c(2,6,10), strand=c('+','+','-'), stacking="dense", fill=c('red','blue', 'black'))
png(file=paste("test","png", sep="."))
plotTracks(track)
dev.off()

```

Note, I am not interested in alternative software recommendations, I would like to solve this with gviz if it is possible. Also asked this here.

# Answer
The primary author of Gviz confirmed to me via email that this is not possible:

> Gviz will only give you the traditional linear tracks. There are a bunch of solutions available for visualizing circular genomes, most of which are more interactive in nature and not implemented in the rather static R graphics context. GGplot2 has some support for circular axes, too. If you want to stay within the R environment I’d start there. Ggbio is another option providing a somewhat higher-level interface.
> 1 votes

---
Tags: r, bioconductor, gviz
---


title: Sequence alignment using Markov Model

text:
# Question
Title: Sequence alignment using Markov Model
I am learning about applying Markov model to sequence alignment. The prof says that the transition probabilities from a gap-residue alignment to a residue-gap alignment and vice versa are both 0. Is there any biological/mathematical reason behind this statement? Why are the (X,Y) and (Y,X) cell 0? This is a lecture slide of lecture 1, week 4 of the "Bioinformatics: Introdcution and Methods" course on coursera.

# Answer
> 6 votes
If I understand your question correctly, then I think for case of pairwise alignment, there is a simple explanation.

I believe the key insight is that: *a mismatch should always score better than a gap*.\*

This follows biologically since the insertion/deletion (indel) rate is roughly 1/10th that of the substitution rate (i.e. the occurrence of single nucleotide changes), at least in vertebrates. (This varies across the tree of life but I think the substitution rate virtually always exceeds the indel rate.)

To understand why this matters, consider an example:

```
ATG-AG
ATGT-G

```

This is an 'impossible alignment' under the probabilities you gave since here we have a transition from a gap-residue alignment to a residue-gap.

However, under our assumption that mismatches are more likely biologically than indels, the correct alignment should be:

```
ATGAG
ATGTG

```

Indeed, the latter does look like a better alignment.

This also follows for more complex examples, so this:

```
ATG--AAG
ATGTT-AG

```

Becomes this:

```
ATG-AAG
ATGTTAG

```

(Or this:

```
ATGA-AG
ATGTTAG

```

)

* Strictly, I mean a substitution should score better than an indel (with the associated gap opening and extension penalties). In fact, for the assumption to always be true, a run of mismatches should still score worse than a single indel. This may not always be a correct assumption, consider this example below, is the true alignment case 1) or 2) or something else? Or is in fact a global alignment bad here and this should be split into 2 local alignments? Is there a likely biological mutational event that could explain this? I ask these questions just to point out it is not black-and-white, I don't have clear answers

1)

```
CGTACGTAGAGGAATGCCCCCCCCC--------AGCAACGTAGCAT
CGTACGTAGAGGAATG---------TTTTTTTTAGCAACGTAGCAT

```

2)

```
CGTACGTAGAGGAATGCCCCCCCCCAGCAACGTAGCAT
CGTACGTAGAGGAATGTTTTTTTT-AGCAACGTAGCAT

```

# Answer
> 1 votes
All Chris\_Rands said is correct: you set the probability of $X\to Y$ and $Y\to X$ to 0 to forbid adjacent insertions/deletions in the alignment. A lot of textbooks including some classical ones use this rule, but in fact, the rule is questionable. It is easier to see this from Smith-Waterman alignment under the affine gap penalty, which is largely the non-probabilistic view of paired HMM.

With the affine gap penalty, a gap of length $k\>0$ is scored as $$ g(k)=-(d+k\cdot e) $$ where $d\ge0$ is the gap open penalty and $e\>0$ is the gap extension penalty. Suppose we are using a simple scoring matrix where a mismatch gets $-b$, $b\>0$. We may see an insertion immediately followed by a deletion (and vice versa) if $b\>2e$. It is actually not so difficult for this to happen. For example, for the human-mouse alignment (see the blastz paper), $e=30$ and $b$ is ranged from 31 to 125. It is possible that an $X\to Y$ transition is preferred in the alignment.

Theoretically speaking, it makes more sense to consider immediate transitions between insertions and deletions. In practice, though, the difference between allowing/disallowing such transitions is probably minor most of time.

EDIT: on Chris' example

```
CGTACGTAGAGGAATGCCCCCCCCC--------AGCAACGTAGCAT
CGTACGTAGAGGAATG---------TTTTTTTTAGCAACGTAGCAT

```

If we use a scoring matrix with $b\>2e$ but disallow adjacent ins-to-del transitions, we will probably end up with an alignment like

```
CGTACGTAGAGGAATGCCCCCCCCC-------AGCAACGTAGCAT
CGTACGTAGAGGAATG--------TTTTTTTTAGCAACGTAGCAT

```

This alignment score will be lower.

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


# Question
Title: AttributeError: module 'scater' has no attribute 'normalize'
I am using `rpy2` to run normalization using `scran` package. After calling `computeSumFactors` I want to do the actual normalization calling `normalize()` function from `scater` package:

```
def normalize(counts):
utils = importr('utils')
utils.chooseCRANmirror(ind=1)

packnames = ('Matrix', 'limSolve')
names_to_install = [x for x in packnames if not isinstalled(x)]
if len(names_to_install) > 0:
    utils.install_packages(StrVector(names_to_install))

base = importr('base')
base.source('http://www.bioconductor.org/biocLite.R')
biocInstaller = importr('BiocInstaller')

biocInstaller.biocLite('SingleCellExperiment', suppressUpdates=True)
biocInstaller.biocLite('scran', suppressUpdates=True)
biocInstaller.biocLite('scater', suppressUpdates=True)

sceLib = importr("SingleCellExperiment")
scranLib = importr("scran")
scaterLib = importr("scater")
named_list = robjects.r.list(counts=counts)
sce = sceLib.SingleCellExperiment(named_list)

try:
    clusters = scranLib.quickCluster(sce)
    sce = scaterLib.computeSumFactors(sce, cluster=clusters, positive=True)
except rpy2.rinterface.RRuntimeError:
    print("fewer cells than the minimum cluster size, so no quickCluster is run")
    sce = scranLib.computeSumFactors(sce, positive=True)

----> sce = scaterLib.normalize(sce, logExprsOffset=1)

```

Everything runs fine except the last line which gives the error:

> AttributeError: module 'scater' has no attribute 'normalize'

This is strange since that is what I am supposed to do actually according to docs.

I suspect that maybe I am loading a wrong `scran` version? I do not know though neither how to check the version number, nor how to install a specific one. Any suggestions would be greatly appreciated.

I tried the following:

```
sce = scaterLib.normalizeSCE(sce, log_exprs_offset=1)

```

And it is working. I am afraid though that it is a wrong way to go.

I tried to get the normalized matrix but it retrieved the initial raw counts one:

```
sce = scaterLib.normalizeSCE(sce, log_exprs_offset=1)
mtx_norm_r = robjects.r.assay(sce)
mtx_norm = np.array(mtx_norm_r)

```

So, something is definitely wrong here. Either I am retrieving the normalized matrix in the wrong way or the normalization was not run correctly. I tried the following method for retrieval but no luck:

```
mtx_norm_r = sceLib.normcounts(sce)

```

> RRuntimeError: Error in assay(object, i = exprs\_values) : 'assay(, i="character", ...)' invalid subscript 'i' 'normcounts' not in names(assays())

I found a way to get the version numbers using `getNamespaceVersion()` function:

```
robjects.r.getNamespaceVersion('scater')

```

The versions are:

```
scater - 1.8.4
scran - 1.8.4
SingleCellExperiment - 1.2.0

```

They are the newest as I understand.

# Answer
I was able to get to the normalized values using `logcounts` function in the following way:

```
mtx_norm_r = sceLib.logcounts(sce)
mtx_norm = np.array(mtx_norm_r)

```

Since I used `log_exprs_offset=1` the field that I needed to use is `logcounts`, otherwise it should be `normcounts` which I found here.
> 1 votes

---
Tags: r, python, scran, rpy2
---


title: K mean clustering issue

text:
# Question
Title: K mean clustering issue
So I have expression data set with 4 healthy sample and 4 disease, to determine the number of K I used mclust which i ran on the list of differentially expressed genes so I got around 7 cluster from Mclust which im using as an input for my K-mean .

Now my data set is like all the differentially expressed genes it contains both UP and down regulated genes so basically it may be UP in my healthy and Down in patient and vice-versa ,4 likely possibilities.

So now I have this doubt that from that I'm getting 7 cluster.Is that result is correct ?

# Answer
There is generally no such thing as correct or incorrect when it comes to clustering. When you specify $K$, the program did it's statistical algorithms. Unless the the program had a bug, your results were statistically valid.

You will need to check the clustered labels visually. Only your biological knowledge can tell you the optimal $K$.
> 1 votes

# Answer
7 clusters may be correct in that it's the optimal number according to the Bayesian information criterion (or whatever value you had mclust optimize), but that doesn't make it necessarily the most sensible value. Have a look at your data and visually look at what might make sense. It's likely that the clusters are subdividing up and down regulated genes into various subgroups, which may or may not be useful, depending on your goals.
> 1 votes

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


# Question
Title: Error in annotating row names as data points in a PCA plot
I made a PCA plot and was trying to plot the eigenvenctors in R so that each data point is actually the sample name. All the samples are in column 1. The following are the R commands I am using:

```
    fn <- "Plot.txt"
    evecdat <- read.table(fn, col.names=c("Sample", "PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8", "PC9", "PC10", "pc11","pc12", "pc13", "pc14", "pc15","pc16", "pc17", "pc18", "pc19", "pc20", "pop", "color"))
    jpeg("Plot1.jpeg", width = 20, height = 20, units = 'in', res = 300)
    identify(evecdat$PC1, evecdat$PC2  = NULL, labels = seq_along(evecdat$Sample), pos = TRUE, n = length(evecdat$Sample), plot = FALSE, atpen = FALSE, offset = 0.5, tolerance = 0.25, …)

    plot(evecdat$PC1, evecdat$PC2, col=as.vector(evecdat$color), xlab="PC1", ylab="PC2", pch="24", cex.lab=2.5, cex=1.5, cex.axis=1)
    legend("topright", ncol=1, legend=c("X", "Y", "Z"), col=c("Blue", "Red", "Green"), pch="24", cex=2, bty = "o")
    while (!is.null(dev.list())) dev.off()

```

Here, all the row names are under the Sample column name that I defined earlier.

# Answer
Maybe try

```
plot(evecdat$PC1, evecdat$PC2, xlab="PC1", ylab="PC2", pch="24", cex.lab=2.5, cex.axis=1, type="n")
text(evecdat$PC1, evecdat$PC2, evecdat$Sample, col=as.vector(evecdat$color))

```

This should plot the sample names as text (taken from the `evecdat$Sample` column), colored as specified in the `evecdat$color` column. See `?text` for details.
> 2 votes

---
Tags: bioconductor, ggplot2, pca
---


# Question
Title: Genomic coordinates of predicted binding pocket
I would like to know the genomic coordinates of a proteins ligand-binding pocket. For example, on this PDB page you can see that there is xray crystallography data for how CYP2C9 interacts with two ligands. I would like to know which amino acids form the space on the protein that binds the ligands, then translate (de-translate?) that back to the nucleotide sequence. Is anyone familiar with either a database that has this information or how to extract the ligand interaction data from the PDB files?

Thank you

# Answer
It turns out there is a tool that does this beautifully. I used the `proteinToGenome` function from the ensembldb package.

The code looks like this:

```
# Define the Ensemble database object for the
edbX <- filter(EnsDb.Hsapiens.v75, filter = ~ seq_name == "X")
# Set the first and last amino acid you want to convert and the ensembl protein id
GAGE10_prt <- IRanges(start = 5, end = 9, names = "ENSP00000385415")
# Perform the conversion.  GAGE10_prt is an object that has the coordinates
GAGE10_prt <- proteinToGenome(GAGE10_prt, edbX)

```
> 3 votes

# Answer
That PDB page includes a publication reference and abstract, containing the following text:

> Mammalian CYP450s recognize and metabolize diverse xenobiotics such as drug molecules, environmental compounds and pollutants.
>
> ...
>
> Here we describe the crystal structure of a human CYP450, CYP2C9, both unliganded and in complex with the anti-coagulant drug warfarin.

Here's the genomic context for CYP2C9, with Ensembl page here.

The binding pocket images from the paper here and here indicate that the following amino acids are potentially involved in the interaction:

* Arg97
* Phe100
* Ala103
* Val113
* Phe114
* Asn217
* Pro367
* Phe476

Using the filtering from the UCSC table browser (genome: Human; group: Genes & Gene Predictions; Basic table; Region: genome; Identifiers: CYP2C9; Output format: CDS fasta alignment), I get the following amino acid sequences respectively for the exons in CYP2C9. I have added in the protein reference range to the end, with the help of the R `cumsum` function applied to the region lengths, `cumsum(c(56,54,50,54,59,47,63,47,61))`:

```
>ENST00000260682.7_hg38_1_9 56 0 0 chr10:94938683-94938850+ [1-56]
MDSLVVLVLCLSCLLLLSLWRQSSGRGKLPPGPTPLPVIGNILQIGIKDISKSLTN
>ENST00000260682.7_hg38_2_9 54 0 1 chr10:94941858-94942020+ [57-110]
LSKVYGPVFTLYFGLKPIVVLHGYEAVKEALIDLGEEFSGRGIFPLAERANRGF
>ENST00000260682.7_hg38_3_9 50 1 1 chr10:94942192-94942341+ [111-160]
GIVFSNGKKWKEIRRFSLMTLRNFGMGKRSIEDRVQEEARCLVEELRKTK
>ENST00000260682.7_hg38_4_9 54 1 0 chr10:94947779-94947939+ [161-214]
ASPCDPTFILGCAPCNVICSIIFHKRFDYKDQQFLNLMEKLNENIKILSSPWIQ
>ENST00000260682.7_hg38_5_9 59 0 0 chr10:94949108-94949284+ [215-273]
ICNNFSPIIDYFPGTHNKLLKNVAFMKSYILEKVKEHQESMDMNNPQDFIDCFLMKMEK
>ENST00000260682.7_hg38_6_9 47 0 1 chr10:94972104-94972245+ [274-320]
EKHNQPSEFTIESLENTAVDLFGAGTETTSTTLRYALLLLLKHPEVT
>ENST00000260682.7_hg38_7_9 63 1 0 chr10:94981183-94981370+ [321-383]
AKVQEEIERVIGRNRSPCMQDRSHMPYTDAVVHEVQRYIDLLPTSLPHAVTCDIKFRNYLIPK
>ENST00000260682.7_hg38_8_9 47 0 1 chr10:94986033-94986174+ [384-430]
GTTILISLTSVLHDNKEFPNPEMFDPHHFLDEGGNFKKSKYFMPFSA
>ENST00000260682.7_hg38_9_9 61 1 0 chr10:94988847-94989028+ [431-491]
GKRICVGEALAGMELFLFLTSILQNFNLKSLVDPKNLDTTPVVNGFASVPPFYQLCFIPVZ

```

And these nucleic acid sequences:

```
>ENST00000260682.7_hg38_1_9 168 0 0 chr10:94938683-94938850+
ATGGATTCTCTTGTGGTCCTTGTGCTCTGTCTCTCATGTTTGCTTCTCCTTTCACTCTGGAGACAGAGCTCTGGGAGAGGAAAACTCCCTCCTGGCCCCACTCCTCTCCCAGTGATTGGAAATATCCTACAGATAGGTATTAAGGACATCAGCAAATCCTTAACCAAT
>ENST00000260682.7_hg38_2_9 163 0 1 chr10:94941858-94942020+
CTCTCAAAGGTCTATGGCCCTGTGTTCACTCTGTATTTTGGCCTGAAACCCATAGTGGTGCTGCATGGATATGAAGCAGTGAAGGAAGCCCTGATTGATCTTGGAGAGGAGTTTTCTGGAAGAGGCATTTTCCCACTGGCTGAAAGAGCTAACAGAGGATTTG
>ENST00000260682.7_hg38_3_9 150 1 1 chr10:94942192-94942341+
GAATTGTTTTCAGCAATGGAAAGAAATGGAAGGAGATCCGGCGTTTCTCCCTCATGACGCTGCGGAATTTTGGGATGGGGAAGAGGAGCATTGAGGACCGTGTTCAAGAGGAAGCCCGCTGCCTTGTGGAGGAGTTGAGAAAAACCAAGG
>ENST00000260682.7_hg38_4_9 161 1 0 chr10:94947779-94947939+
CCTCACCCTGTGATCCCACTTTCATCCTGGGCTGTGCTCCCTGCAATGTGATCTGCTCCATTATTTTCCATAAACGTTTTGATTATAAAGATCAGCAATTTCTTAACTTAATGGAAAAGTTGAATGAAAACATCAAGATTTTGAGCAGCCCCTGGATCCAG
>ENST00000260682.7_hg38_5_9 177 0 0 chr10:94949108-94949284+
ATCTGCAATAATTTTTCTCCTATCATTGATTACTTCCCGGGAACTCACAACAAATTACTTAAAAACGTTGCTTTTATGAAAAGTTATATTTTGGAAAAAGTAAAAGAACACCAAGAATCAATGGACATGAACAACCCTCAGGACTTTATTGATTGCTTCCTGATGAAAATGGAGAAG
>ENST00000260682.7_hg38_6_9 142 0 1 chr10:94972104-94972245+
GAAAAGCACAACCAACCATCTGAATTTACTATTGAAAGCTTGGAAAACACTGCAGTTGACTTGTTTGGAGCTGGGACAGAGACGACAAGCACAACCCTGAGATATGCTCTCCTTCTCCTGCTGAAGCACCCAGAGGTCACAG
>ENST00000260682.7_hg38_7_9 188 1 0 chr10:94981183-94981370+
CTAAAGTCCAGGAAGAGATTGAACGTGTGATTGGCAGAAACCGGAGCCCCTGCATGCAAGACAGGAGCCACATGCCCTACACAGATGCTGTGGTGCACGAGGTCCAGAGATACATTGACCTTCTCCCCACCAGCCTGCCCCATGCAGTGACCTGTGACATTAAATTCAGAAACTATCTCATTCCCAAG
>ENST00000260682.7_hg38_8_9 142 0 1 chr10:94986033-94986174+
GGCACAACCATATTAATTTCCCTGACTTCTGTGCTACATGACAACAAAGAATTTCCCAACCCAGAGATGTTTGACCCTCATCACTTTCTGGATGAAGGTGGCAATTTTAAGAAAAGTAAATACTTCATGCCTTTCTCAGCAG
>ENST00000260682.7_hg38_9_9 182 1 0 chr10:94988847-94989028+
GAAAACGGATTTGTGTGGGAGAAGCCCTGGCCGGCATGGAGCTGTTTTTATTCCTGACCTCCATTTTACAGAACTTTAACCTGAAATCTCTGGTTGACCCAAAGAACCTTGACACCACTCCAGTTGTCAATGGATTTGCCTCTGTGCCGCCCTTCTACCAGCTGTGCTTCATTCCTGTCTGA

```

Associating the pocket protein amino acid numbers in exons 2,3,5,7,9 with genomic sequence using these FASTA sequences is left as an exercise for the reader.
> 1 votes

---
Tags: proteins, genomics, drugs
---


title: using SNPs to identify mixed samples

text:
# Question
Title: using SNPs to identify mixed samples
Is there a way to identify mixed samples based on SNPs?

Example input (table of genotypes for multiple samples):

```
|      | sample1 | sample2 | sample3 |
|------|---------|---------|---------|
| rs1  | AA      | BB      | AB      |
| rs2  | AA      | BB      | AB      |
| rs3  | AA      | BB      | AB      |
| rs4  | AA      | BB      | AB      |

```

Desired output: Sample relationships/distances. In this example, sample1 and sample2 are distinct, but sample3 is likely a mix of two (biological or technical).

Based on this simple example, the calculation seems doable, but maybe it's impossible in more complex realistic scenarios.

There is a nice tool peddy that does something conceptually similar, but it expects pedigree information as well as genotype frequencies for both the input samples and population.

# Answer
Yes. The general term for identifying subsets from a mixed group is deconvolution.

This is easier to do when there is haplotype-level information available (e.g. sequenced reads that span multiple SNPs), but it is possible to get some indication of sample mixes from unlinked SNPs alone.

Deconvolution is helped a lot by having proportional data in the input (e.g. frequencies of genotypes, rather than a single genotype). Most SNPchip platforms have multiple probes per variant, and can provide some indication of variant proportion when digging down into the raw data. Deconvolution is also a lot easier with some knowledge about the component mixes. It's much easier to simulate mixes than to guess at possible mixes.

Here's a paper about the application of deconvolution in forensics.

Here's a paper about a program called CIBERSORT which carries out a deconvolution of immune populations based on gene counts:

Here's another paper which discusses the idea in more general terms.

You might find more information on one of these Biostars posts:
> 2 votes

---
Tags: snp, genotyping
---


# Question
Title: Extract certain fields using from GenBank file using Bash script
I need to extract certain lines from a GenBank file on Linux.

Input file:

```
LOCUS       XM_012540068            1713 bp    mRNA    linear   MAM 11-JAN-2018
DEFINITION  PREDICTED: Sarcophilus harrisii cytochrome c oxidase assembly
            protein COX15 homolog (LOC100925550), mRNA.
ACCESSION   XM_012540068
VERSION     XM_012540068.2
DBLINK      BioProject: PRJNA167725
KEYWORDS    RefSeq; corrected model.
SOURCE      Sarcophilus harrisii (Tasmanian devil)
  ORGANISM  Sarcophilus harrisii
            Eukaryota; Metazoa; Chordata; Craniata; Vertebrata; Euteleostomi;
            Mammalia; Metatheria; Dasyuromorphia; Dasyuridae; Sarcophilus.
COMMENT     MODEL REFSEQ:  This record is predicted by automated computational
            analysis. This record is derived from a genomic sequence
            (NW_003816462.1) and transcript sequence (GEDN01126617.1) annotated
            using gene prediction method: Gnomon, supported by EST evidence.
            Also see:
                Documentation of NCBI's Annotation Process

            On Jan 11, 2018 this sequence version replaced XM_012540068.1.

            ##Genome-Annotation-Data-START##
            Annotation Provider         :: NCBI
            Annotation Status           :: Full annotation
            Annotation Name             :: Sarcophilus harrisii Annotation
                                           Release 102
            Annotation Version          :: 102
            Annotation Pipeline         :: NCBI eukaryotic genome annotation
                                           pipeline
            Annotation Software Version :: 8.0
            Annotation Method           :: Best-placed RefSeq; Gnomon
            Features Annotated          :: Gene; mRNA; CDS; ncRNA
            ##Genome-Annotation-Data-END##

            ##RefSeq-Attributes-START##
            assembly gap :: added 136 transcript bases to patch genome assembly
                            gap
            ##RefSeq-Attributes-END##
PRIMARY     REFSEQ_SPAN         PRIMARY_IDENTIFIER PRIMARY_SPAN        COMP
            1-136               GEDN01126617.1     1-136
            137-327             AEFK01003589.1     1071-1261
            328-450             AEFK01003589.1     2007-2129
            451-637             AEFK01003589.1     2869-3055
            638-805             AEFK01003589.1     6071-6238
            806-887             AEFK01003591.1     516-597
            888-1042            AEFK01003591.1     1431-1585
            1043-1156           AEFK01003591.1     4257-4370
            1157-1713           AEFK01003591.1     6191-6747
FEATURES             Location/Qualifiers
     source          1..1713
                     /organism="Sarcophilus harrisii"
                     /mol_type="mRNA"
                     /db_xref="taxon:9305"
                     /chromosome="1"
                     /map="unlocalized"
                     /sex="female"
     gene            1..1713
                     /gene="LOC100925550"
                     /note="The sequence of the model RefSeq transcript was
                     modified relative to its source genomic sequence to
                     represent the inferred CDS: added 136 bases not found in
                     genome assembly; Derived by automated computational
                     analysis using gene prediction method: Gnomon. Supporting
                     evidence includes similarity to: 10 ESTs, 11 Proteins, and
                     100% coverage of the annotated genomic feature by RNAseq
                     alignments, including 5 samples with support for all
                     annotated introns"
                     /db_xref="GeneID:100925550"
     CDS             53..1288
                     /gene="LOC100925550"
                     /note="The sequence of the model RefSeq protein was
                     modified relative to its source genomic sequence to
                     represent the inferred CDS: added 84 bases not found in
                     genome assembly"
                     /codon_start=1
                     /product="cytochrome c oxidase assembly protein COX15
                     homolog"
                     /protein_id="XP_012395522.2"
                     /db_xref="GeneID:100925550"
                     /translation="MMRLLFPPLKAVLGSPGVRLQAAGGAHRCGCHYGQKLSLLQRQH
                     STLSEVALQSGRGAVSVPTKAAQRIMGRWLLVCSGTVAGAVILGGVTRLTESGLSMVD
                     WHLIKGMKPPTSQEEWEAEFQKYQQFPEFKILNHNMTLAEFKFIWYMEYSHRMWGRIV
                     GLAYVLPTVYFWRKGWLTHGMKGRVLALCGLVCFQGLLGWYMVKSGLEEKPDSYDIPR
                     VSQYRLAAHLGSALVLYCASLWTGLSLLLPRDKLPETRKLLQLRRFAHGTAGLVFLTA
                     LSGAFVAGLDAGLVYNSFPKMGDHWIPEDLFTFSPVLRNVFENPTMVQFDHRVLGIAT
                     VTAITGLYFFSRRIPLPRRAKMAVTTLLAVAYAQVALGISTLLTYVPTPLAAVHQSGS
                     LALLSIALWLMNELRRIPK"
ORIGIN
        1 gtggggcgac cgagggaggc cgctgcgagg ctgtgacttt gagactcgcg cgatgatgcg
       61 attgctgttt cctcctctga aagcggtgct gggcagtccg ggcgtccggc tgcaggctgc
      121 tgggggagct caccggtgcg gttgccatta tgggcaaaag ctctctctct tgcaaagaca
      181 gcacagtact ctgtctgaag tagctttgca gtctggaagg ggagcagtgt ctgtccccac
      241 aaaggctgcc cagcgcatta tgggtcgttg gctcttggtc tgcagcggaa cagtagcagg
      301 agcggtgatt cttggtggag tgactaggtt gacagagtct gggctctcca tggtagattg
      361 gcacttgata aaaggaatga aaccacccac aagccaagaa gaatgggaag cagagttcca
      421 aaaataccag cagtttccag aatttaaaat cttgaaccat aacatgacat tagcagaatt
      481 caagtttatc tggtacatgg aatactcaca cagaatgtgg ggtcggattg tgggcctagc
      541 atatgtcctt cccactgtct acttttggag aaagggctgg ctcacccatg gcatgaaagg
      601 tcgtgttctt gccctctgtg gtttagtctg tttccaggga ctgttggggt ggtatatggt
      661 gaagagtgga ttagaagaaa agccagattc ttatgacatc cctcgggtga gccagtatcg
      721 tctggccgct cacttgggct ctgcccttgt cctctactgt gccagcttgt ggactggact
      781 ctccctcttg cttcctcgag acaagttgcc agaaacccga aagctcttgc agctgaggcg
      841 atttgctcat gggacagcag gcctggtttt ccttacagct ctttcaggag cctttgtggc
      901 aggcctggat gctggacttg tatacaactc ctttcctaag atgggagacc actggatccc
      961 agaagatctc tttaccttct ctccagttct tagaaatgtc tttgagaatc caaccatggt
     1021 gcagtttgat catcgagtcc tgggaatcgc tacagtcaca gccattactg gactttactt
     1081 cttctccagg agaatccctc ttcctcgaag ggccaaaatg gcagtgacaa cattgctggc
     1141 tgtggcctat gcacaggtgg ccttgggcat tagtactcta ctaacgtatg ttccaactcc
     1201 tttggctgca gttcaccagt cagggtccct cgctttgctc agcattgccc tttggcttat
     1261 gaacgaactg aggcgaatcc caaagtaatg cccaggcagg tatttttctg gactcttccc
     1321 aggcttttga gagctctgta gagggcatga ggcttccctc ctgttgaatt aaagaattga
     1381 tgtggtagac agatcagaag atcaagaatg aagttactga tctgaggatt ctttttcaga
     1441 ctaaacaaac tactcttctt ccataagttg tcacttgctt tttaatgtgc tattaagatc
     1501 ccttccaact ttcactactc tgattctttt gatacagcaa taccactgat aggtctgtgt
     1561 ctcagagttg ggagggaggg aaggaaaaaa aaaacttttt gtacaaacgt acaaaactat
     1621 agtagttcct tttgtgggaa ctaaaataag gaaattttct ttctttttta tttattcgtt
     1681 ttcaacattt acttttataa gatttagatt ttc
//

```

My code:

```
tr "\"" " " < XM_012540068.2.gb | egrep -v "/db_xref= taxon" | awk -F= '/GeneID:|(gene|product|protein_id|translation)=/{print $2}' | sort -u

```

What I get:

```
 GeneID:100925550
 LOC100925550
 MMRLLFPPLKAVLGSPGVRLQAAGGAHRCGCHYGQKLSLLQRQH
 XP_012395522.2
 cytochrome c oxidase assembly protein COX15

```

What I want:

```
 GeneID:100925550
 LOC100925550
 MMRLLFPPLKAVLGSPGVRLQAAGGAHRCGCHYGQKLSLLQRQH...(sequence goes on...)
 XP_012395522.2
 cytochrome c oxidase assembly protein COX15 homolog

```

Preferably, it has to be a mastercode for every GenBank-file. For a project I've got to loop through 78 GenBank-files, and extract the same data. I asked a similar question earlier, but the solution there doesn't work when the protein's name spans multiple lines.

# Answer
Here is a script for achieving your result not so good but giving the desired output

```
#!/bin/bash
egrep -v "/db_xref= taxon"  $1| awk -F= '/GeneID:|(gene|protein_id)=/{print $2}' | sort -u > file1
sed -n '/product/,/protein_id/p' $1|grep -v protein_id |sed 's/.*=//' >file2
awk '/translation/,/ORIGIN/{if(/translation|ORIGIN/) next; print}' $1 >file3
head -3 file1 && cat file2 && printf "\"" && cat file3 |tr -d " \t\n\r" &&   echo && rm -rf file1 file2 file3

```

Save it in a file names as script1 and executing it as

```
`./script5 filename`

```

Gives following output

```
"GeneID:100925550"
"LOC100925550"
"XP_012395522.2"
"cytochrome c oxidase assembly protein COX15
                 homolog"
"STLSEVALQSGRGAVSVPTKAAQRIMGRWLLVCSGTVAGAVILGGVTRLTESGLSMVDWHLIKGMKPPTSQEEWEAEFQKYQQFPEFKILNHNMTLAEFKFIWYMEYSHRMWGRIVGLAYVLPTVYFWRKGWLTHGMKGRVLALCGLVCFQGLLGWYMVKSGLEEKPDSYDIPRVSQYRLAAHLGSALVLYCASLWTGLSLLLPRDKLPETRKLLQLRRFAHGTAGLVFLTALSGAFVAGLDAGLVYNSFPKMGDHWIPEDLFTFSPVLRNVFENPTMVQFDHRVLGIATVTAITGLYFFSRRIPLPRRAKMAVTTLLAVAYAQVALGISTLLTYVPTPLAAVHQSGSLALLSIALWLMNELRRIPK"

```
> 0 votes

---
Tags: proteins, genbank, bash
---


# Question
Title: What does "ALT+ cancer cells", the plus sign mean?
In many scientific papers, there is a plus sign:

> ... T-circles have also been detected in ALT+ (telomerase-negative) cells9 and ...

What does the "+" mean? Does it mean differentially expressed?

# Answer
> 5 votes
Here are a couple of papers ( and ) that explain ALT.

> Telomere maintenance is a hallmark of cancer as it provides cancer cells with cellular immortality. A significant fraction of tumors uses the alternative lengthening of telomeres (ALT) pathway to elongate their telomeres and to gain an unlimited proliferation potential. Since the ALT pathway is unique to cancer cells, it represents a potentially valuable, currently unexploited target for anti-cancer therapies.

The ALT+ denotes cell lines that are positive for the ALT pathway.

---
Tags: terminology
---


title: Fastqc- Per Base Sequence Quality

text:
# Question
Title: Fastqc- Per Base Sequence Quality
I am trying to figure out how to interpret the "Per Base Sequence Quality"? What does Position in reads (bp) mean? Also in order to draw this box-plot graphics, more than one quality scores are needed so we can sort them and calculate the mean, median and quartiles. In this diagram, how these values were calculated for each position?

# Answer
> 2 votes
There is a pretty good description here. The position in the reads means literally the position in the reads, so position 2 means the quality of all bases on position 2 in the fastq file (all reads). In the description link, you can find how the boxplots and blue and red lines are calculated (and what they mean). They are all based on the total quality values of all bases on a certain position, all bases means all bases of all the reads on a particular position.

---
Tags: ngs, sequence-analysis, fastq
---


title: How to get the product of a CDS

text:
# Question
Title: How to get the product of a CDS
I need the name of the protein in /product="protein\_name" using bash commands. Beware, there is a lot of whitespace between lines.

```
FEATURES             Location/Qualifiers
     source          1..1266
                     /organism="Sarcophilus harrisii"
                     /mol_type="mRNA"
                     /isolate="Salem"
                     /db_xref="taxon:9305"
                     /chromosome="1"
                     /map="unlocalized"
                     /sex="female"
                     /country="Australia: Taronga Zoo, Sydney"
     gene            1..1266
                     /gene="UROS"
                     /note="Derived by automated computational analysis using
                     gene prediction method: Gnomon. Supporting evidence
                     includes similarity to: 10 Proteins, and 100% coverage of
                     the annotated genomic feature by RNAseq alignments,
                     including 5 samples with support for all annotated
                     introns"
                     /db_xref="GeneID:100933602"
     CDS             60..908
                     /gene="UROS"
                     /codon_start=1
                     /product="uroporphyrinogen-III synthase isoform X1"
                     /protein_id="XP_003755055.2"
                     /db_xref="GeneID:100933602"
                     /translation="MRIPTSFSRGAWWKMKVLLLKDAKDDDSGLDPYLQELGSCGFEA
                     TLIPVLAFEFLSLPSFSEKLSHPEAYGGLVFTSPRAVEALQLSLERAGQLEVWKKSLK
                     EKWRLKPAYVVGQATASLVKKVGLAAEGEKCGNGEKLAEYICSRETPATLPLLFPCGT
                     LKGEVLPKTLKQKGILLESLTVYQKVQHPRMQDSLKSYFSTQGVPAGVTFFSPSGLTY
                     SLPHIQELSGASFDRIKFAAIGPTTARALASAGIPVSCTAEKPCPQDLAAGLQKALQP
                     SGPSSS"

```

Right now I have this:

```
while read p;
do cat "${p}.gb" | tr "\"" " "| egrep -v "/db_xref= taxon"| egrep "/gene=|/db_xref=|/product=|/protein_id=" | sort | uniq | awk '{print($2)}END{print("")}'
done <XM_codes.txt

```

Which results in this:

```
GeneID:100933602
UROS
uroporphyrinogen-111
XP_003755055.2

```

My output should be:

```
GeneID:100933602
UROS
uroporphyrinogen-III synthase isoform X1
XP_003755055.2

```

# Answer
If you have GNU `grep` (default on Linux), you can do this with a single `grep` command using Perl Compatible Regular Expressions (PCREs):

```
$ grep -oP 'GeneID:\d+|/(gene|product|protein_id)="\K[^"]+' file | sort -u
GeneID:100933602
uroporphyrinogen-III synthase isoform X1
UROS
XP_003755055.2

```

The `-o` flag tells `grep` to only print the matching portion of the line and the `-P` enables PCREs which let us use `d+` for "one or more digits" and `\K` which means "discard whatever was matched before the `\K`".

Note that I'm not using `cat` here. The command from your question is what is known as a UUoC (useless use of cat) and is a surefire way of annoying any \*nix geeks around ;).

---

If you don't have GNU `grep`, you can do something like this:

```
$ awk -F= '/GeneID:|(gene|product|protein_id)=/{print $2}' file |
    tr -d '"' | sort -u
GeneID:100933602
uroporphyrinogen-III synthase isoform X1
UROS
XP_003755055.2

```
> 0 votes

# Answer
If you change the delimiter to something besides a space in the `tr` command (`tr "\"" "\t"`) then you can use that same delimiter in the `awk` command (`-F "\t"`). In my example below I used a tab, also I filtered the taxon prior to the `tr` so that you don't have to worry about that staying in sync if you change the delimiter.

```
 egrep -v "/db_xref=\"taxon"| tr "\"" "\t" | egrep "/gene=|/db_xref=|/product=|/protein_id=" | sort | uniq | awk -F "\t" '{print($2)}END{print("")}';

```

results:

```
GeneID:100933602
UROS
uroporphyrinogen-III synthase isoform X1
XP_003755055.2

```
> 3 votes

# Answer
The solutions here are more robust than mine since I rely on the structure of the file, but if you want a solution that is really easy to understand:

```
$grep -A 5 "CDS" file | grep -v "CDS" | grep -v "codon_start" | cut -f 2 -d \"
UROS
uroporphyrinogen-III synthase isoform X1
XP_003755055.2
GeneID:100933602

```

I take CDS line and 5 lines bellow, then I kick out CDS and codon\_start (lines we are not interested in) and then I just bring what is between quote using cut. I also report the values in the order as they appear in the original file.
> 0 votes

---
Tags: proteins, linux, genbank, bash
---


title: Boxplot of kmean cluster

text:
# Question
Title: Boxplot of kmean cluster
```
 gene   H4       H3       H2       H1       L8       L7       L6       L5       L4       L3       L2       L1          clust
AAGAB  7.540712 7.296325 7.030218 7.408974 8.269541 7.758752 9.253013 8.611002 8.705045 8.940146 8.855133 8.488776     3
ABCA1  7.112951 6.213973 7.200973 6.784828 8.044299 9.832347 8.558342 7.134570 8.993938 5.794208 8.271960 9.179501     5
ABCA6  1.822567 2.209167 3.379618 2.813433 3.180290 3.741516 3.060194 2.078485 6.536346 6.881403 5.927874 5.278854     2

```

This is my output after doing kmean clustering so at the end I have my cluster number assigned to each gene. Now I want to show how the expression of my sample varies across each cluster, I am not sure how to transform the data and do it in `ggplot2` as of now i doing something like this in base R plot

```
boxplot(da$H4 ~ k$cluster,
        xlab='Cluster', ylab='Genes',
        main=' Genes by Cluster')

```

For each sample I am checking, is there any way to show that with all the samples together, any suggestion or help would be highly appreciated.

# Answer
> 1 votes
Assuming `da` is a data-frame with the cluster as factor you can do the following:

```
ggplot(da) +
  geom_boxplot(aes(cluster, H4)) +
  xlab("Cluster") +
  ylab("Genes") +
  ggtitle("Genes by Cluster")

```

---
Tags: r, visualization
---


# Question
Title: Aligning sequence and comparing it against primer
I am looking to show how a primer is consistent among some genomic data. I have a primer of about 23bp and looking to compare it to about 5000 genomic sequences of 10kb. I am unsure how to format it the way I need to, since I am just trying to show the nucleotide differences of the primer to the genetic sequences. Here is what I need:

```
>     1. Cut out the area that my primer is located and about 20bp down each end.
>     2. Show only the bases that are different from my primer in my analysis.
>      ex: Primer: -----------ATGTGGAAGCAAATATCAAATGA---------
>          Gene:   ATGACCATACG----C--------------T---ATCGTAGGG
>                  ATGAGCATACC-----A----T--------T---TTCGAACGC

```

The data I am using is all dengue sequences (all serotypes) and the primer with the following code: ATGTGGAAGCAAATATCAAATGA.

I am using Biostrings, msa, and seqinr. I use ncbi to get the genetic sequences and using FASTA files.

Thanks!

Edit: Solution: I used MAFFT to align the data, working very efficiently and put it into a fasta file. Then I cut the nucleotides I did not need by: 1. Finding the part i wanted and cutting 10 nucleotides before it and after it. 2. Putting only the \[start:end\] into a new list, which I will then manipulate.

However, I am still having difficulty finding a way of showing my data nicely because I can't find a suitable package. Any advice?

Edit: Hi all, Just wanted to show the code I made to solve the problem.

```
    #So the we will be using a dataframe that has the sequences as a string
    #Then I created a temporary dataframe to store each sequence, with the primer at the top

    tempDF <- data.frame(matrix(primer1, nrow = 1, ncol = 1),
    stringsAsFactors = FALSE)

    for(i in sequenceAsString){
      #Temporary vector for each sequence after compared to primer.

      temporary1 <- ""
      temporary1 <- paste(temporary1, substring(i,1,11), sep = "")

      for(j in 12:(nchar(i)-20)){
        if(substring(i,j,j) == substring(toString(reversePrimer1), j-12,j-12)){
          temporary1 <- paste(temporary1, "-", sep = "")
        }else{
          temporary1 <- paste(temporary1, substring(i,j,j), sep = "")
        }
      }
      temporary1 <- paste(temporary1, substring(i,nchar(i)-19,nchar(i)), sep = "")
      tempDF <- rbind(tempDF,temporary1)
    }

```

So, I essentially take the sequence I chose apart and compare each character to that of the reverseprimer. I then replace the ones that are the same with an '-'

Thanks everyone for their help!

# Answer
> 3 votes
If I understand you correctly, after your most recent edit, you have 5,000 sequences that are 63 bp long and in each you simply want to replace the nucleotides from 21 through 43 with `-` unless that nucleotide does not match the primer. Since `-` is normally used to represent gaps, I recommend using `.` instead for matching nucleotides.

I know you're looking for an R solution, but I don't usually use R for string manipulations and don't have time at the moment to figure it out. So, assuming this might be urgent for you, below is a Python solution. This uses your example data which shows 11 nucleotides before the primer and 9 after it (so you'll need to adjust it to work with your actual case in which 20 nucleotides are before and after the primer). In Python, strings are arrays of characters, by the way.

Here are your example sequences are in the file `seqs.txt`:

```
ATGACCATACGATGTCGAAGCAAATATCAATTGAATCGTAGGG
ATGAGCATACCATGTGAAAGCTAATATCAATTGATTCGAACGC

```

Here is the code that gives you the expected output:

```
#!/bin/env python
import sys
primer = 'ATGTGGAAGCAAATATCAAATGA'

print('-' * 11 + primer + '-' * 9 )
with open('seqs.txt','r') as f:
    for gene in f:
        sys.stdout.write(gene[0:11])
        for gene_index in range(11,34):
            primer_index = gene_index - 11
            if gene[gene_index] == primer[primer_index]:
                sys.stdout.write('.')
            else:
                sys.stdout.write(gene[gene_index])
        print(gene[34:])

```

Running this gives you:

```
-----------ATGTGGAAGCAAATATCAAATGA---------
ATGACCATACG....C..............T...ATCGTAGGG
ATGAGCATACC.....A....T........T...TTCGAACGC

```

---
Tags: r, bioconductor, sequence-alignment, primer
---


title: Convert SRA to FastA

text:
# Question
Title: Convert SRA to FastA
I'm trying to get the FastA files for some accessions (like NC\_001416.1). I did not managed to find an FTP server or direct link to these files (I want to get it from command line with `wget`, not from a web browser). But I found an "equivalent" file in ftp://ftp.ncbi.nlm.nih.gov/sra/refseq/NC\_001416.1

This file is an SRA file that should be processed with the SRA-toolkit. In particular, fastq-dump seems to be the tools of choice.

I ran

```
fastq-dump --fasta NC_001416.1

```

and I got a `NC_001416.1.fasta` file. But the file contains the sequence I'm looking for in chunks of 5000 bases. The IDs in the output file are:

```
>NC_001416.1.1 length=5000
>NC_001416.1.2 length=5000
>NC_001416.1.3 length=5000
>NC_001416.1.4 length=5000
>NC_001416.1.5 length=5000
>NC_001416.1.6 length=5000
>NC_001416.1.7 length=5000
>NC_001416.1.8 length=5000
>NC_001416.1.9 length=5000
>NC_001416.1.10 length=3502

```

I need that sequence under a single ID:

```
>NC_001416.1 length=48502

```

I cannot find the proper parameters to supply to `fastq-dump` to achieve the desired result.

Does anyone know how could I get the FastA file I'm expecting? How to process properly the SRA file or where I could download the FastA file directly.

# Answer
Finally, I found an alternative to the SRA translation: a link that works! For those of you interested in knowing how to download FastA files from NCBI using an accession number, try the following link:

```
https://www.ncbi.nlm.nih.gov/search/api/sequence/${accession}/?report=fasta

```

Using `wget` to download the accession used as example:

```
wget -O NC_001416.1.fasta "https://www.ncbi.nlm.nih.gov/search/api/sequence/NC_001416.1/?report=fasta"

```
> 5 votes

# Answer
Here is a Biopython solution via `Entrez` utilies, simply copy the code into a python script and execute it:

```
from Bio import Entrez, SeqIO

Entrez.email = 'me@email.com'

def get_fasta_from_ids(ids_list):
    handle = Entrez.efetch(db='nucleotide', id=','.join(ids_list), rettype='fasta')
    for record in SeqIO.parse(handle, 'fasta'):
        yield record

SeqIO.write(get_fasta_from_ids(['NC_001416.1']), 'NC_001416.1.fa', 'fasta')

```

You can add more ids to the input list to retrieve multiple records at once.
> 2 votes

# Answer
You can either do this directly as you show in your answer or, for a more sophisticated and flexible approach, use NCBI's `edirect` tool:

```
esearch -db nucleotide -query 'NC_001416.1' | efetch -format fasta > NC_001416.1.fa

```
> 1 votes

# Answer
You can add the `--origfmt` option, and see what comes out. Either you will have what you want, or perhaps they uploaded the sequences that way.

See: https://ncbi.github.io/sra-tools/fastq-dump.html
> 0 votes

---
Tags: fasta, format-conversion, data-download, sratoolkit
---


# Question
Title: DESeq2 complicated design - effect of replicated samples
I have RNAseq data from a relatively complicated experimental design with variables = genotype, treatment, time, and batch. I have 2 biological replicates for each genotype/condition, however unfortunately in the first iteration of the experiment one sample had poor RNA quality and is unusable. To get around this lack of a replicate, we repeated the experiment for the two biological replicates in the condition of the lost sample, and did a new library prep/sequencing from this repeated experiment. The repeated experiment obviously lead to a batch effect, and although it isn't the main driver of variance, it is a significant factor. I would like to model or remove the batch effect from the experiment using DESeq2 if possible, but after reading this similar post: can I model technical replicates in DESeq2? I am concerned about inflating the power of the test for those conditions that were replicated.

**Here is a diagram of the experimental design to make it more easy to understand:**

My "colData" for the DESeq2 object looks like this:

My previous thought was to simply collapse the genotype+time+treatment variables into a single combined factor "condition" and then model the batch effect using a design of "~batch + condition". However this doesn't really address the concern of inflating the power of the tests due to technical replication.

Should I Just collapse the counts of these replicated conditions even though the batch affect will affect the combined count values? Or should I use a different tool than DESeq2 ? Or is there a way for me to remove the batch effect prior to collapsing the counts for replicates?

**What would you guys recommend?**

# Answer
> 2 votes
I agree with @devon-ryan that the design matrix looks reasonable, i.e. that it's okay to treat the new samples as biological replicates of the original batch, and that the batch number should be included in that matrix (assuming DESeq2 doesn't complain about it) to filter out most of the batch effects.

DESeq2 tends to be conservative in its differential expression calculations. While it's true that more replicates in one sample will increase the number of differentially-expressed genes that are statistically significant (because the biological variation can be better modelled), I wouldn't expect this to change the overall picture.

However, whole-genome differential expression testing should not be used on its own to demonstrate biological relevance. I recommend as a first-pass sanity check that people compare any differential expression results to what is seen on a genome browser to identify any quirky things going on (e.g. expression across only a part of the gene). Ideally, results should be confirmed with targeted analyses on the same (or similar) biological samples, either at the transcript level (e.g. qPCR) or at the protein level (e.g. Western Blot). When aggregating evidence from multiple analyses, the specific way in which one of those analyses was done is less important.

Aside: It'd obviously be better \[in the future\] to make sure that your sample replicates are both robust to drop-out, and allow for outliers to be noticed, reducing the batch-effect problem. In practise this means a *minimum* of four replicates per informative condition, but because funders / researchers are well known for only choosing the minimum, I'd just state that six replicates per sample as what is needed to produce good results.

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


title: Using Broad Institute RNAseq pipeline

text:
# Question
Title: Using Broad Institute RNAseq pipeline
There is a Github report for Broad's pipeline on RNA-Seq. It has a docker build file and WDL scripts.

However, I don't see how any of the WDL script is being used in the docker (Dockerfile). It looks like Broad's docker runs on the standard command-line tools, but not via the WDL language. How do the Broad bioinformaticians take advantage of the WDL scripts?

# Answer
You need to follow the readme of this folder to use the pipeline, there you will find instructions about how to align, mark duplicates, quantificate, and aggregate the results.

If you want to execute the WDL scripts you can follow this guide using Cromwell. Basically use: `java -jar Cromwell.jar run myWorkflow.wdl --inputs myWorkflow_inputs.json`
> 2 votes

---
Tags: rna-seq, docker, wdl
---


# Question
Title: Understanding SNP coding for association analysis
I'm working on a project about detecting SNP association with a disease. As I understand, SNP is a single variation of the nucleotide that occurs for more than 1% of the population. So, if a gene is AT**T**G and there's a variation AT**A**G then the SNP is **A**. However, I couldn't connect this idea with the dataset in hand. The rows in my dataset represent each patient and the columns contain SNP information. For example:

```
ID  exm355 exm615
1     T_T   A_C
2     T_T   C_C
3     A_T   C_C

```

I have no idea why the SNP columns contain 2 nucleotides(T\_T, A\_T, A\_C, C\_C). As the definition of SNP, I thought it should show only the variant nucleotide or am I misunderstanding anything? How could I interpret T\_T or C\_C and how could I know which nucleotide is the variation from the common ones in the population?

The source says that the case and control subjects were genotyped by the Australian Genome Research Facility using the Illumina HumanExome BeadChip-12v1\_A, which exonic content consists of $\>$ 250 000 markers representing diverse populations and a range of common conditions. An individual was duplicated to test genotyping reliability and quality. The identity by descent (IBD) matrix between all pairs of individuals was estimated after linkage disequilibrium (LD) pruning and used for quality control and for the mixed linear models analyses.

# Answer
Humans are diploid. This means they have two copies of each locus (with the exception of sex chromosomes in males). Your two alleles for each individual represent the two copies of this variant locus.
> 3 votes

---
Tags: snp
---


# Question
Title: RNA velocity: competing explanations for variable ratio of spliced to unspliced transcript
In "RNA velocity of single cells", La Manno et al. look at ratios between spliced and unspliced mRNA as a way of estimating the velocity of a cell through transcriptome space.

In checking for a universal steady-state spliced-unspliced ratio, they write "We found that the steady-state behaviour of most genes across a wide range of cell types was consistent with a single fixed \[ratio between spliced and unspliced\] (Extended Data Fig. 3a–c). However, 11% of genes showed distinct slopes in different subsets of tissues (Extended Data Fig. 3d, e), suggesting tissue-specific alternative splicing (Extended Data Fig. 3f) or degradation rates."

Tissue-specific alternative splicing or degradation rates certainly could account for these patterns, but how do La Manno et al. rule out that some tissues are simply not in a steady state? The figure they cite (screenshotted below) uses data from the Tabula Muris adult mouse cell encyclopedia. In an adult mouse, not all cells would not be expected to be in steady state -- for example, T-cell progenitors from the thymus and epithelial cells in many organs undergo continual replenishment.

# Answer
If cells aren't in a steady state then one wouldn't expect to see the results in (e), where there are two slopes. One would instead expect to see results more like Plekha3 in panel (c), where there's more variability around the the steady state line (or perhaps no steady state line at all).
> 1 votes

# Answer
In addition to Devon Ryan's data-driven answer, I've found a clue in the text. La Manno et al have not used the full Tabula Muris: the figure legend indicates that they used only terminally differentiated cells.
> 1 votes

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


# Question
Title: How to find GEO data sets using Drug Bank ID and bioDBnet
I want to find all experiments in GEO that are associated with a drug (for example tolvaptan). Is there any quick and scalable way to to this? I want to query more than 100 drugs.

I tried to use bioDBnet to map Drug Bank ID to look up data sets, but I don't know what output format should I select.

# Answer
You could take a look at the `rentrez` package for R. You can use the `entrez_search` function to get the GEO IDs for the studies related to your query, and then the `entrez_summary` function for the results. So, for example:

```
x <- rentrez::entrez_search(db="gds", term=paste0("tolvaptan AND gse[ETYP]"),
                            use_history = T)
res <- rentrez::entrez_summary(db="gds", id=x$ids)
str(res)

```

Where `res` will be a list containing objects of the `esummary` class. Then, you can use the following line to get the GSE accession numbers related to tolvaptan:

```
lapply(res, `[[`, "gse")

```

So you could put this function in a loop replacing "tolvaptan" with your drugs:

```
drugs <- c("drugA", "drugB", "drugC")
for(drug in drugs){
    x <- rentrez::entrez_search(db="gds", term=paste0(drug, " AND gse[ETYP]"))
    res <- rentrez::entrez_summary(db="gds", id=x$ids)
    lapply(res, `[[`, "gse")
}

```

Take a look at the `rentrez` documentation for more details.
> 3 votes

# Answer
One of the problems I've found with GEO and even ArrayExpress is that it seems like there are tons of false positives that come up in the search. Additionally, there may be a ton of results that you miss because you didn't include the synonym of the drug.

In order to make sure I have all the information relating to name, etc, you can write a script which downloads the result of a search in PubChem (this can be done with drugbank as well I think) and then parse the search result, extracting important information such as synonyms of drug name.

Finally, take those synonyms and use them in an Entrez (EUtils) search for the experiments. Then try to filter the results to keep only what you're actually looking for. One bad way to do this is just by making sure the experiment summary contains the term your looking for. Another way to potentially do this, for example if there is a large experiment which tests serveral different factors is to programmatically search through associated files to ensure they contain a desired keyword.

You'll get different experiment types and so the files associated with the results will depend on the type of experiment in the search result.

When you get have a list of IDs of the experiments you want, you can feed them through to GEOquery in Bioconductor.
> 2 votes

# Answer
I would recommend using GEOmetadb from Bioconductor. GEOmetadb is a simple wrapper that will download an up-to-date SQLite database with the metadata for all entities in GEO. You can then use your favorite SQL client to query the database quickly, since it's offline, refining your search as you go.
> 0 votes

---
Tags: database, gene, drugs
---


# Question
Title: How can the Autocovariances, autocorrelations, and autocorrelation coefficients be calculated from a Protein Amino Acid Sequence?
Given a normal protein sequence with the 20 standard amino acids, how can the 'Autocovariances', 'autocorrelations', and 'autocorrelation coefficients' of the sequence be calculated?

What is meant by these terms in the context of a protein sequence?

What are the appropriate algorithms and formulas?

Is there any existing package/function to calculate these values?

# Answer
Protein autocorrelation is similar to usual autocorrelation of time series (that is, function $r(k)$ measuring correlation between values lagged by $k$), except that instead of time dimension we slide along the backbone of amino acid chain. The "dependent variable" can be any quantifiable property of the amino acids - then the autocorrelation measures similarity in that particular property.

See e.g. Fernandez et al., 2008, for more extensive background: https://www.tandfonline.com/doi/full/10.1080/08927020802301920?scroll=top&needAccess=true
> 2 votes

---
Tags: proteins, statistics, features
---


# Question
Title: GATK documentation for required depth to reliably call heterozygous mutation in diploid organism?
I'm looking for official GATK documentation (or a recent manuscript) that defines a general recommendation/requirement for sequencing depth to reliably call a heterozygous point mutation in a diploid organism (WGS). In this case, I'm working with humans. Specifically, suppose we're looking at position X in the genome and we want to classify (True/False) whether there was enough coverage to detect a heterozygous point mutation.

What coverage would that be? I'm curious about both: (1) UnifiedGenotyper and (2) HaplotypeCaller

I know general practice is to have 30x coverage genomewide, but I don't see anything in GATK's documentation (i.e., Best Practices) about required coverage.

Really appreciate your help.

Note: I also submitted a question to the GATK team end of last week.

# Answer
> 2 votes
**UPDATE:** As an update, Sarah Walker (co-author on the poster) responded to my question on the `GATK` forum. She clarified with the following statement:

> We believe the sites around 30X (and above 150X) are being filtered due to low mapping quality (since it is whole genome so there are many areas that are hard to map), which explains the low sensitivities in these areas. As the error bars show, there are very few variants in these low mapping quality areas.

She then posted the following plot that falls much more inline with what we would all expect. This is a plot showing sensitivity for `NA12878` compared to NIST truth.

---

shlee, from the GATK forums, pointed me to a poster from AGBT 2018 where the authors compared the HiSeqX and NovaSeq instruments. Figure 3 shows sensitivity for SNPs and short INDELs across multiple allelic fractions (AF) from spike-ins. ***I was surprised to find that SNPs where AF = 0.4 only had ~50% sensitivity at ~25x coverage, and only ~75% sensitivity at ~50x coverage***.

The poster is targeted at somatic variation, but data from figure 3 should behave like a standard diploid germline variant because you often see a wide range of allelic fractions from germline variants. They used spike-ins for figure 3. It's not clear to me whether they used HaplotypeCaller in figure 3, but I don't see how they could have used Mutect2, like they did in figures 1 & 2.

Based on this information, **it appears the standard 30x genome-wide coverage is generally insufficient to confidently call many heterozygous mutations, and that ~90-100x is ideal.**

**Am I missing something? Are we really missing that many heterozygous mutations in standard WGS/WES?**

I'll leave this question unanswered for now to allow others to chime in on interpretation, etc. I'd appreciate additional thoughts.

---
Tags: gatk, wgs, coverage
---


title: Seurat with normalized count matrix?

text:
# Question
Title: Seurat with normalized count matrix?
I know that in `Seurat` we have the function `CreateSeuratObject` from which the analysis starts, but it accepts raw count matrix according to the documentation. I have only the already normalized count matrix, so is there a way to work with `Seurat` using normalized data? I would try maybe to feed in the normalized data instead of raw to `CreateSeuratObject`, but I do not know what `Seurat` is doing behind the scenes, whether it is a good idea. Besides, `Seurat` provides by default only one `log` \- normalization method, but I may want to normalize the data by myself with various methods and only then start the analysis with `Seurat` \- that is the other reason why I want to find a way to start from normalized data.

# Answer
> 6 votes
This was addressed by the Seurat developers here:

> if you have TPM counts, I suggest you don't use Seurat::NormalizeData(), since TPM counts are already normalized for sequencing depth and transcript/gene length. Note that Seurat::NormalizeData() normalizes the data for sequencing depth, and then transforms it to log space. If you have TPM data, you can simply manually log transform the gene expression matrix in the object@data slot before scaling the data.

And here:

> If you want to run your own normalization function, please place the data in log scale before placing it into object@data

And here:

> If you're supplying already normalized values to CreateSeuratObject, you should skip the normalization step. You will likely see a warning when you run ScaleData (which is checking whether you normalized using Seurat). You can ignore this or turn it off by setting check.for.norm to FALSE.

This is a lot of quotes, but I found them all to be helpful in different ways.

---
Tags: r, seurat
---


# Question
Title: What is the correct way to map Hi-C data with bwa mem?
## Library Prep

* I have a Hi-C library prepped using an enzyme that cuts at `GATC`, so it leaves `GATCGATC` as the junction sequences. This library was sequenced on a 2x150 PE Illumina run.

## Data Pre-processing

* The reads were adapter and quality trimming with trimmomatic paired-end mode, keeping reads that overlap.
* Finally, the forward and reverse reads were trimmed such that all bases 3' of the junction sequence (`GATCGATC`). For example, the read:
    + `5' AAAATT**GATCGATC**AATTATGGATAA 3'` becomes: `5' AAAATT**GATC 3'`
* In this process, sequences that don't contain the junction sequences were left alone. Trimming junction sequences in this way makes it such that each read pair maps to only two places in the genome. No three-way or four-way chimeric mappings are possible since everything `3'` of the junction sequence was removed.

## My Question

* What is the correct way to now map these trimmed Hi-C reads to a genome using `bwa mem`?

## Progress so far/things I have tried

* In the past I have just tried mapping using `bwa mem` with no special parameters using this command: `bwa mem ref.fa hi-c_f.fq.gz hi-c_r.fq.gz`. One concern is that the computations that `bwa mem` performs on insert sizes in this mode will cause some mappings to be faulty or that read pairs with chromosome-length information will be removed from the alignment because they fall on the tail of the insert size distribution (non-Gaussian, but log decay).
* Another option to consider is mapping each file, `hi-c_f.fq.gz` and `hi-c_r.fq.gz`, independently then adding mate information with `samtools fixmate`. However, I do not know how to do this.
* A last option that I have considered is to use `bwa mem` in paired-end mode, but to use a special flag. In the past I believe that `-5` was the correct flag, but in this post Heng Li discusses that the `-5` flag behavior has changed, as well as a new `-q` flag. In this post Heng Li discusses that he thinks `bwa mem` is still a good option for mapping Hi-C data but does not discuss the proper flags.

## Limitations

* The answer should provide the correct way to map Hi-C using `bwa mem`, but other software suggestions will also be considered.

# Answer
See here, in particular slide #10 from the tutorial:

```
bwa mem -5SPM ref.fa read1.fq read2.fq > out.sam

```

Here, `-SP` disables pairing. The Aiden lab and this paper also use a similar command line. Beware that there are many Hi-C pipelines, but not many are using bwa-mem.
> 4 votes

# Answer
There are some interesting suggestions from the HiCExplorer developers:

> Mates have to be mapped individually to avoid mapper specific heuristics designed for standard paired-end libraries.

```
# map the reads, each mate individually using
# for example bwa
#
# bwa mem mapping options:
#       -A INT        score for a sequence match, which scales options -TdBOELU unless overridden [1]
#       -B INT        penalty for a mismatch [4]
#       -O INT[,INT]  gap open penalties for deletions and insertions [6,6]
#       -E INT[,INT]  gap extension penalty; a gap of size k cost '{-O} + {-E}*k' [1,1] # this is set very high to avoid gaps
#                                  at restriction sites. Setting the gap extension penalty high, produces better results as
#                                  the sequences left and right of a restriction site are mapped independently.
#       -L INT[,INT]  penalty for 5'- and 3'-end clipping [5,5] # this is set to no penalty.

$ bwa mem -A1 -B4  -E50 -L0  index_path \
    -U mate_R1.fastq.gz 2>>mate_R1.log | samtools view -Shb - > mate_R1.bam

$ bwa mem -A1 -B4  -E50 -L0  index_path \
    -U mate_R2.fastq.gz 2>>mate_R2.log | samtools view -Shb - > mate_R2.bam

# build matrix from independently mated read pairs
# the restriction sequence GATC is recognized by the DpnII restriction enzyme

$ hicBuildMatrix --samFiles mate_R1.bam mate_R2.bam \
                 --binSize 10000 \
                 --restrictionSequence GATC \
                 --threads 4
                 --inputBufferSize 100000
                 --outBam hic.bam \
                 -o hic_matrix.h5
                 --QCfolder ./hicQC

```
> 1 votes

---
Tags: read-mapping, bwa, illumina, hi-c
---


# Question
Title: TagReadWithGene missing when using latest version of Drop-seq_tools
I am using Drop-seq\_tools to analyze scRNA-seq in a similar way of this paper. I am using the same data (GSE97930) to get used to Drop-seq\_tools.

On the GEO repo the authors provide the relevant code (GSE97930\_README.txt). However, one tool used in the script is missing: **TagReadWithGene**.

The latest version of Drop-seq\_tools does not have such method (as it should be provided by this library), and there is no version history for the whole package, so I cannot download previous versions.

The best matches for the missing method are **TagReadWithGeneExon** and **TagReadWithInterval**. In the methods of the cited paper they say:

> Reads that mapped to intronic or exonic regions of genes as per the GENCODE gene annotation were recorded.

I think that in this case the tool to use should be **TagReadWithInterval** since the other one should consider only exons.

I cannot find any other documentation about the tools. Do you know/already encountered this problem?

# Answer
> 1 votes
As per this response from Drop-seq\_tools team, at the moment intronic regions cannot be annotated and **TagReadWithGeneExon** should be used instead.

However, in the version of the public repository:

> TagReadWithGeneExonFunction is the old code, for people who want the more simple tagging methodology and rely on it for 3rd party tools, and TagReadWithGeneFunction, which for your purposes replaces TagReadWithGeneExon. \[...\] You’ll see additional options to extract expression for CODING, UTR, INTRONIC reads \[...\]

---
Tags: scrnaseq, gene
---


# Question
Title: How do population genetics people define a population?
How do population genetics people define a population?

Do they define it as a layman will do? say Africans, Americans, or so? Or is there a more scientific way of doing so? For example, I think defining one's population as one's allele frequencies makes more sense, but has this point been discussed previously?

# Answer
It really depends on the study and the ethnicity. Mostly they use self-defined ethnicity and often they will require that three or four grandparents are of that ethnicity. Here's some detailed info on how the ethnicity was defined for the 1000 Genomes populations. You'll see that some populations (such as ASW) require all four grandparents to identify as that ethnicity, but some (such as ACB) only need three, and the cultural sensitivities in Japan meant that JPT participants weren't questioned directly but told that the study was only interested in people with Japanese grandparents. The mixed American super-population (AMR) is an exception to this rule due to the extreme admixture of ethnicity in these regions, and is based more in geography than ethnicity – this means that the individuals had to be from the place and have grandparents from the place.
> 2 votes

# Answer
A population is defined as a group of people (or organisms of the same species) living together in a geographical region. However, for human genetic studies, it is wise to use ethnicity as well (so group of people with the same ethnic background). See here how they define population in the 1000 genomes project.
> 3 votes

---
Tags: phylogenetics, genetics
---


# Question
Title: Check if a set of genes fall more on a specific cluster than other clusters
I have a network which I clustered into 11 smaller clusters. I have a set of genes and I would like to check if statistically significant number of these genes fall on any specific cluster. One method which I thought of was to find which percentage of a total cluster were made up of my genes of interest. For example, say 60% of the genes in a cluster were the genes in my list. However, I was hoping that there was a better way to do this. This study on honey bees has done something of the kind.

I would appreciate if there was a way to do something similar. I do not know how to replicate what they have done in the study. Any help is appreciated.

# Answer
> 1 votes
A hypergeometric test is related to the fisher test, in a table with one dimension being some factors and the other one some other factors it is tested if they are independent. In R you con do this:

```
data <- matrix(c(0, 1, 2, 5), ncol = 2, nro = 2)
fisher.test(data)

```

Or you can use the specific function `phyper`, to see the relation between them you can read this question and answer

---
Tags: networks, clustering
---


# Question
Title: How to detect mismatch before mapping in RNA-Seq data
In the methods of this paper, the authors say:

> Reads were then filtered based on quality score in the UMI region. Any read with \>1 low-quality base (phred \<=10) were discarded. **Reads with more than one mismatch in any of the three 8 nt cell barcodes were also discarded**. The cDNA reads (Read 1) were then mapped to either a combined mm10- hg19 genome or the mm10 genome using STAR (44).

I understand how to filter based on mismatch *after* mapping but in this case it is done *before* mapping. Has this to be interpreted as "Reads with *low quality in more than one base* in any of the three 8 nt cell barcodes were also discarded"?

# Answer
> 2 votes
At each step in the SLiTing and pooling the barcodes added are not random, but are taken from a predefined list of possible barcodes.

When the reads arrive, first you filter by reads that have \>1 low quality bases in the UMI region (which is truely random).

Then you compare the sequences in the 3 cell-barcode regions against the list of allowed barcodes for each of those positions. If a read is not within an edit distance of one from one of those barcodes, it is discarded.

Finally the reads that pass these filters can be aligned.

# Answer
> 1 votes
It should not be interpreted as "Reads with low quality in more than one base in any of the three 8 nt cell barcodes were also discarded". They only care if the cell barcodes match what are valid cell barcodes within 1 mismatch, the individual qualities are ignored. This sort of filtering turns out to be pretty common in practice, since you have to write a maximum edit distance into whatever you use to demultiplex by cells (or labeling each read, as was likely done in this case).

# Answer
> 0 votes
Cell barcodes =/= sequence reads. The default settings on Illumina's bcl2fastq work exactly the same way; they allow a read to count for a sample if there is a single mismatch in the index, but not more.

---
Tags: scrnaseq, read-mapping, reads, barcode
---


# Question
Title: What is and how to detect a dephased read
In the methods of this paper, the authors say:

> To simplify analysis, we first removed any dephased reads in our library (last 6 bases of read did not match the expected sequence).

I have read this post and this one, and I think that a dephased read is one (according to the first paper) where the last 6 bases are not detected with good quality. In particular, I think this may be an example of a dephased read:

```
@...
CACACACACCAGTAAAAAANAGTTGCNNCNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
+
AAAAAEEAEEEEEEEEEE6#EAEEEE##E#####################################

```

where the quality of the read drops in the middle. Is that correct? If yes, should they be all of bad quality or at least one? And, if not, how can I detect a dephased read?

# Answer
> 4 votes
Phasing in Illumina data means that the read you get is shifted one base before or after what the sequence really should be. They must be expecting some constant 6 bases at the end of each read, and they are throwing out all reads that don't have those 6 bases. This would include N's, but a small deletion or insertion would also cause the last 6 bases to not match. This is in theory detectable and fixable; the dropSeq protocol for instance catches instances where improperly short cell barcodes shift the UMIs, but these authors look like they decided not to try and bother, and just throw out everything that didn't match their exact expectations.

---
Tags: scrnaseq, reads, splitseq
---


# Question
Title: Calculating bit score: How do you find lambda and K?
To calculate bitscore from score you can use this equation: $S' = (lambda\*S - ln(K)) / ln(2)$

If I am trying to manually calculate the bitscore of an HSP of a pairwise blastn alignment, and I know the alignment score, how do I calculate the parameters K and lambda?

Is there a more simple equation for pairwise nucleotide alignments?

# Answer
> 3 votes
NCBI BLAST uses tables of pre-computed values for lambda/K for gapped protein and nucleotide comparison and their values depend on the amino acid substitution matrix or the substitution scores + gap penalties.

You can find these tables (from line 589 for blastn) in the file `c++/src/algo/blast/core/blast_stat.c` from the blast source package.

---
Tags: blast, bitscore
---


# Question
Title: The visualisation of a list of genes on URD object
Seurat R package has some functions like **FeaturePlot**, **DimPlot** and **DoHeatmap** by which we can plot the expression of a list of genes on cell clusters. I am working with **URD** that likely does not have such options or I can not find that. URD has plotDot but would be clear only for a few number of genes (please look at the picture). Could somebody please help me in plotting the expression pattern of a list of genes in URD? As I am not good in R at all to write my own function or adopt Seurat functions to URD

For example Seurat DoHeatmap gives such a picture by which we can see these set of genes are being expressed on which group of cells

```
DoHeatmap(object ,genes.use=a vector of genes)

```

Or we can plot a set of genes on seurat feature lot to they are being expressed more on which cluster

```
FeaturePlot(object = pbmc_small, features.plot = 'PC1')

```

In URD only by

```
plotDot(object.6s.mnn, genes = a vector of genes, clustering="Infomap-60")

```

We can plot a set of genes but as you are being in the first picture that would not be clear much

My question is how I can plot a vector of genes on URD object

for example Seurat DimPlot function can visualise the expression pattern of a group of genes on cell clusters but when I am trying to do the same with URD, returns this error

```
> plotDim(axial.6somite, c("DDB_G0291149","DDB_G0291029"), plot.title="SHHA (Notochord marker)")
Error in names(data) <- cells.use :
  invalid to use names()<- on an S4 object of class 'dgCMatrix'
In addition: Warning messages:
1: In if (label.type == "meta" | (label.type == "search" & label %in%  :
  the condition has length > 1 and only the first element will be used
2: In if (label.type == "group" | (label.type == "search" & label %in%  :
  the condition has length > 1 and only the first element will be used
3: In if (label.type == "sig" | (label.type == "search" & label %in%  :
  the condition has length > 1 and only the first element will be used
4: In if (label.type == "gene" | (label.type == "search" & label %in%  :
  the condition has length > 1 and only the first element will be used
>
>

```

Actually only one gene can be visualised by this function on URD

# Answer
I just added URD clustering labels for each cel to seurat object and then I did DoHeatmap by using group.by="URD clustering labels"
> 1 votes

---
Tags: r, scrnaseq, seurat
---


title: Spearman correlation between two genes

text:
# Question
Title: Spearman correlation between two genes
I have FPKM data for three genes like below:

```
     Samples      gene1 gene2 gene3
TCGA-2Y-A9GS-01A   0.9   7.4  1.0
TCGA-2Y-A9GT-01A   0.8   1.0  0.3
TCGA-2Y-A9GU-01A   0.6   2.0  0.2
TCGA-2Y-A9GV-01A   1.2   0.5  0.1
TCGA-2Y-A9GW-01A   3.8   2.1  0.4
TCGA-2Y-A9GX-01A   2.3   2.0  1.5

```

I used `cor.test`

```
cor.test(~ gene1 + gene2, data = df2, method="spearman", continuity=FALSE, conf.level=0.95)

    Spearman's rank correlation rho

data:  gene1 and gene2
S = 5686100, p-value = 6.083e-09
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho
0.2984045

```

I have a warning message which I didn't see before.

```
Warning message:
In cor.test.default(x = c(0.9, 0.8, 0.6, 1.2, 3.8, 2.3, 3.8, 0.4,  :
  Cannot compute exact p-value with ties

```

Do I need to care about this warning message? Is it good using FPKM data for correlation?

For plotting I used `ggscatter`.

```
ggscatter(data, x = "gene1", y = "gene2",
          add = "reg.line", conf.int = TRUE,
          cor.coef = TRUE, cor.method = "spearman",
          xlab = "gene1", ylab = "gene2")

```

Is this fine or I need to use any `log` for the scatter plot?

# Answer
Since Spearman is a rank-based test, it relies on you being able to accurately decide on the ranking of your observations by some metric (usually the magnitude of the numbers). If two observations have identical values (are *tied*), then they cannot be definitively ranked. Since the ranks are not unique, exact p-values cannot be determined.

e.g. in your data, for gene 2:

```
                  gene2 rank
TCGA-2Y-A9GS-01A    7.4    1
TCGA-2Y-A9GT-01A    1.0    5
TCGA-2Y-A9GU-01A    2.0   3=
TCGA-2Y-A9GV-01A    0.5    6
TCGA-2Y-A9GW-01A    2.1    2
TCGA-2Y-A9GX-01A    2.0   3=

```

By way of example, it's possible to do the test with either `TCGA-2Y-A9GU-01A` or `TCGA-2Y-A9GX-01A` ranked number 3 - meaning a definitive p-value cannot be determined.

If there are only a few ties, you probably don't need to worry too much. If your FPKMs (or better, TPMs) were calculated with greater precision, then ties would be less likely. Is there a reason you only have these numbers to 1 decimal place?

Log or linear scale doesn't make a huge difference (and neither is 'wrong'). If you choose to plot log counts, you need to consider what happens if you have zeroes in your data.
> 4 votes

---
Tags: rna-seq, fpkm, correlation
---


# Question
Title: How to iterate protein sequences using amino acids?
I was working on certain program using python and I have been using a protein aligned sequence file in two formats, phylip (.phy) and clustal (.aln). Example clustal file:

```
CLUSTAL 2.1 multiple sequence alignment

    Homo            -------------------------------------------MEMQDLTSPHSRLSGSS
    Macaca          -------------------------------------------MEMQDLTSPHSRLSGSS
    Canis           REWSSARPERSKGRRKPVDAAAVSAVQTSQTSSDVAVSSSCRSMEMQDLTSPHSRLSGSS
    Mus             -------------------------------------------MEMQDLTSPHSRLSGSS
    Rattus          -------------------------------------------MEMQDLTSPHSRLSGSS
                                                               *****************

```

Example phylip file:

```
     3   120
Apis       MPWIKCLGGG GGKARRGKPL SVSQPIHLLV KSEFEPQCHV FRTKQLRKPR
Branchiost ---------- ---------- ---------- ---------- ----------
Nematostel ---------- ---------- ---------- ---------- ----------

           PCHLCHQAVI KQASCCRVCK YICHKTCEDK GARFRDTSQQ IHKQWQADPA
           ---------- ---------- ---------- ---------- ----------
           ---------- ---------- ---------- ---------- ----------

```

I want to navigate through the file in a way that ;

1. If the columns of the file have same amino acid, it saves and prints the amino acid that exists in the column.
2. I need to know how can I search the columns by using the amino acids, if particular group of amino acids exist in that column or not. For example, if 'STA'or 'HY' or 'FVLIM' or 'NDEQHK' (with no '-') exists in any of the columns, then put a '@' in the file at the end of that column. Following is the code I have been trying to manipulate the file with the following program:

Code below

```
alignment = open ("example.aln")
a = b = c = 0
for col in (alignment):
    num = len(set(col))
        if num == 1 and col[0] != '-':
            print (num)
        a += 1
        elif num > 1 and '-' not in col:
        print (num)
        b += 1 elif '-' in col:  # assumes 1 or more dashes
        c += 1
        print('Number of columns with the same amino acid: {}\n'
        'Number of columns with at least 2 amino acids (no gaps): {}\n'
        'Number of columns with one/more gaps: {}'
        .format(a, b, c))

```

These variables return the occurrences of amino acids in the columns of the file but I don't understand how to search using the amino acids. What is the best way to iterate over the file and write into a new file with results using python? Should I use python switch statement?

# Answer
In response to question 1, save the following as a python3 script (e.g., `aln_parser.py`):

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

from Bio import AlignIO
import sys

aln = sys.argv[1]

algnmnt = AlignIO.read(open(aln, 'r'), 'clustal')

for col in range(0, algnmnt.get_alignment_length()):
    proteins = set(algnmnt[:,col])
    if len(proteins) == 1:
        print(f"Postion {col}: protein {''.join(proteins)}")

```

Run the script to query your alignment (`clustal.aln`).

`python3 aln_parser.py clustal.aln`

To save the output to file:

`python3 aln_parser.py clustal.aln > results.txt`

Requires `python3` and `Biopython` (`pip3 install biopython`).
> 3 votes

---
Tags: python
---


# Question
Title: How can I convert SCF trace files to ABI files?
I have a large collection of sequence trace files which are in SCF file format (binary files that in addition to the string of DNA bases also contain the electropherogram of the sample and quality information about the base calls).

I have an existing BioPython pipeline that accepts files in AB1 format (pdf), so I need to convert my SCF files to AB1 files.

I found a question on Biostars that asks for the opposite: converting from SCF to AB1.

The accepted answer there was to use `convert_trace` from the Staden package.

```
./convert_trace -out_format scf < trace.ab1 > trace.scf

```

So I tried the opposite:

```
./convert_trace -out_format abi < trace.scf > trace.ab1

```

This did indeed produce a file with `*.ab1` extension, but upon importing it to my BioPython pipeline, I got the error

```
File should start ABIF, not '\xaeZTR'

```

So apparently, `convert_trace` converted my files to ZTR format, even altough I specified `-out_format abi`.

Closer reading of the manual states that `[convert_exe] can write CTF, EXP, PLN, SCF and ZTR formats.` (so indeed no ABI), but it confuses me why `ABI` is still a valid output format.

So, my question is, how do I convert SCF files into ABI files?

# Answer
> 4 votes
Apparently `convert_trace` does not do a good parameter checking and silently sets ZTR as default if it does not recognise a valid output format.

I don't know of any freely available command-line converters in that specific direction. Even at the height of Sanger sequencing (late 90s, early 2000s), pipelines usually tried to reduce complexity and footprint by converting AB1 to anything else, often SCF (ZTR came too late and never got enough traction).

That being said, I dug out an SCF I had from 1997, loaded it into Geneious, and re-exported it as .ab1. I was a bit astonished to see that this worked. So Geneious (or similar packages from other companies) may be your best bet for a quick, one-off thing.

Depending on your needs, you might want to think about other ways: either expand BioPython to read SCFs, or maybe to expand the Staden io\_lib to be able to write .ab1 (James still maintains that package and is extremely helpful).

Good luck.

# Answer
> 1 votes
DNA Baser Assembler will convert ABI to SCF (and other formats). Just double click the ABI file you want to convert then go to the "Samples" menu and choose "Save current samples as...". Plus, it can automatically trim the low quality ends.

There is another tool in that package what will convert all ABI files at batch: Chromatogram Explorer. There is a Lite version that is free, but the conversion to SCF might only work for the first 100 files (if I remember correctly).

Related to the comment posted by BioGeek about the Genious license: DNA Baser Assembler does not require a license key. It will work right away.

---
Tags: format-conversion
---