# Question
Title: How to obtain amino acid context at a given position
I'm trying to find the amino acid context around a particular position in a protein (5 aa amino acid window on each side).

I used to have some code that works and then biomaRt stopped working. I think some of these issues would be fixed in later versions of biomaRt, but I would have to install a new version of R, which would require me to install a new version of my OS.

```
require('biomaRt')
require('Biostrings')
hs_ensembl_mart =
useMart(biomart="ENSEMBL_MART_ENSEMBL",dataset="hsapiens_gene_ensembl")

gene_names=c('KRAS')
aa_affected_loc=12

biomart_query_results=getBM(attributes=c('hgnc_symbol', 'entrezgene','coding'),
  filters = 'hgnc_symbol',
  values = gene_names,
  mart = hs_ensembl_mart)

preceding_peptide= substr(as.character(translate(DNAString(biomart_query_results[1,'coding']))),aa_affected_loc-5,aa_affected_loc-1)

trailing_peptide= substr(as.character(translate(DNAString(biomart_query_results[1,'coding']))),aa_affected_loc+1,aa_affected_loc+5)

```

Does anyone have a workaround for this?

# Answer
> 2 votes
It turned out that I needed to reinstall OSX to upgrade my R version and biomaRt install.

# Answer
> 2 votes
You could also save your "biomart\_query\_result" object after the first time it is generated. Then you can re-use it without having to depend on the biomaRt API functionality.

# Answer
> 1 votes
Would this be useful to you?

```

   library(EnsDb.Hsapiens.v86)
   library(ensembldb)
   mygene <- genes(EnsDb.Hsapiens.v86, filter= GenenameFilter("KRAS"))

```

# to get the protein info and protein sequence for each transcript

proteins(EnsDb.Hsapiens.v86, filter= GenenameFilter("KRAS"))

Also have a look at functions proteinToTranscript proteinToGenome

```
> x<- proteins(EnsDb.Hsapiens.v86, filter= GenenameFilter("KRAS"), return.type = "AAStringSet")
> x
  A AAStringSet instance of length 4
    width seq                                    names
[1]   188 MTEYKLVVVGAGGVGKSA...KDGKKKKKKSKTKCVIM ENSP00000308495
[2]    75 MTEYKLVVVGAGGVGKSA...KDGKKKKKKSKTKCVIM ENSP00000452512
[3]   189 MTEYKLVVVGAGGVGKSA...KEEKTPGCVKIKKCIIM ENSP00000256078
[4]    43 MTEYKLVVVGAGGVGKSA...HFVDEYDPTIEVLLSSM ENSP00000451856

```

There are also functions exons, exonsBy, transcriptsBy, cdsBy, fiveUTRsByTranscript and threeUTRsByTranscript. These functions return a GenomicRanges object and inherit all of the methods in that package too.

```
x<-  genes(EnsDb.Hsapiens.v86, filter= GenenameFilter("KRAS"))
start(x)
end(x)
strand(x)
length(x)
flank(x)

```

You can manipulate the sequence selected as an IRanges object

```
ranges(x)   # gives an IRanges object
narrow(ranges(x), start=10, end=20)  # select the 10th position

```

```
> ranges(x)
IRanges object with 1 range and 0 metadata columns:
                      start       end     width
                  <integer> <integer> <integer>
  ENSG00000133703  25204789  25250936     46148
> narrow(ranges(x), start=10, end=19)
IRanges object with 1 range and 0 metadata columns:
                      start       end     width
                  <integer> <integer> <integer>
  ENSG00000133703  25204798  25204807        10

```

It is worth having a look at these vignettes. Let me know if you've questions or post to the Bioconductor support forum ;-)

---
Tags: r, biomart
---


# Question
Title: multiple COSMIC id for the same mutation
I would like to know why there are different entries for the same mutation in COSMIC, for example COSM6954941 and COSM12833 both refer to ERBB4 c.908C\>A.

In this specific case the field **Gene name** in the two entries is different, COSM6954941 report ERBB4\_ENST00000402597 while COSM12833 report ERBB4, but there are many cases in which the **Gene name** is the same, for example COSM1361 and COSM1219.

# Answer
> 3 votes
Harry Boutselakis addressed the question via the customer support of the COSMIC web site.

---

> Dear Ivan, Thank you for using COSMIC.
>
> > I would like to know why there are different entries for the same mutation in COSMIC, for example COSM6954941 and COSM12833 both refer to ERBB4 c.908C\>A.
> >
> > In this specific case the field Gene name in the two entries is different, COSM6954941 report ERBB4\_ENST00000402597 while COSM12833 report ERBB4,
>
> The mutations are annotated onto two different isoforms of ERBB4. So this annotation refers to two different transcripts. In these cases the consequence of the mutation on the mRNA and polypeptide is the same reflected by the syntax.
>
> > but there are many cases in which the Gene name is the same, for example COSM1361 and COSM1219.
>
> This case the consequence one mutation on the mRNA is unknown , but the amino acid change is the same on both cases.
>
> The information stored in COSMIC is curated from experimental data, depending on the evidence and experimental methodology used, the mutations are differentiated accordingly. I hope this answers your question.
>
> Kind regards,
>
> Harry

---
Tags: variation
---


# Question
Title: How can I install bioconductor-gviz and use it in jupyter notebook?
I tried to install gviz in a conda environment, but that library seems to be incompatible to python and r. I tried to setup a clean environment using

```
conda create -n r -c conda-forge r-essentials jupyter

```

and then add the library with:

```
source activate r
conda install -c bioconda bioconductor-gviz

```

getting

> UnsatisfiableError: The following specifications were found to be in conflict: \- atk

I remove atk, but now I get:

> UnsatisfiableError: The following specifications were found to be in conflict: \- bioconductor-gviz \- r-bindr

Does anyone manage to use gviz from within an jupyter notebook?

I also tried to install gviz from a running R-notebook with:

```
install.package('gviz')

```

> Warning in install.packages : package ‘gviz’ is not available (for R version 3.4.3)

Same when I try 'bioconductor-gviz'.

# Answer
> 3 votes
Did you try following the installation instructions?

Try this in R:

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

```

# Answer
> 4 votes
```
# Install the BiocManager package first
install.packages("BiocManager")
library("BiocManager")

# Then install the Gviz package
BiocManager::install("Gviz")

```

# Answer
> 1 votes
Use BiocManager::install()

```
# Code
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("Gviz", version = "3.8")

```

However, If R \< 3.5.0 , please use the following for installing Bioconductor packages.

```
BiocInstaller::biocLite("Gviz")

```

---
Tags: bioconductor, gviz, bioconda, jupyter
---


title: R package installation in r studio

text:
# Question
Title: R package installation in r studio
I have written a package for my personal use and my other colleagues. I have uploaded the files and source package to different websites. I am able to install the source package in my computer. But my colleagues cannot install it either by github.com (https://github.com/furqan915/Epi-gene) or by sourceforge.net (https://sourceforge.net/projects/epi-gene/). Common error is package cannot be installed and it does not add to user library. I am unable to find out the problem. Kindly help me out.

# Answer
Can you use `devtools::install_github()` giving your github username/reposname

If not, you need to check that your R package is passing `R CMD check` or `devtools::check()`. See http://r-pkgs.had.co.nz/check.html and http://kbroman.org/pkg\_primer/pages/check.html for info
> 3 votes

# Answer
A quick workaround is to install the package directly from source:

```
install.packages("~/Download/myPackage.tar.gz", repos = NULL, type="source")

```

You should be able to pass also `.zip` and URLs (other than CRAN).
> 0 votes

---
Tags: r
---


# Question
Title: Extracting genes differentially expressed by Wilcox test
I have 2 cell types (NOF and CAF) cultures individually and also have been co-cultured with tumour like this picture.

https://ibb.co/x8Ty1Wz

If this is the log cpm of these 4 samples

```
> head(onco[,1:4])
             A10       A11        A12        A2
A2M    11.755100 12.145143 12.9787977 12.407221
ABCB11  2.955457  2.278363  0.1229266  1.493507
ABCC2   6.711643  6.974025  0.1229266  7.628936
ABCC6   6.443439  6.094670  6.1389496  7.483394
ABCF1   9.995098 10.045613  9.6627424 10.314842
ABCG2   7.125483  7.448127  6.8669702  8.492895
>

```

I have used log of count per million (cpm) for Kruskal-Wallis rank sum test that says these samples are significantly difference like below;

```
> NOF1=c(NOF1$`onco$A10`)
> NOF2=c(NOF2$`onco$A11`)
> CAF1=c(CAF1$`onco$A12`)
> CAF2=c(CAF2$`onco$A2`)
> my_data <- data.frame(
+     group = rep(c("NOF1", "NOF2","CAF1","CAF2"), each = 719),
+     cell_type = c(NOF1,NOF2,  CAF1,CAF2)
+ )

> head((my_data))
  group cell_type
1  NOF1 11.755100
2  NOF1  2.955457
3  NOF1  6.711643
4  NOF1  6.443439
5  NOF1  9.995098
6  NOF1  7.125483
>

> kruskal.test(cell_type ~ group, data = my_data)

    Kruskal-Wallis rank sum test

data:  cell_type by group
Kruskal-Wallis chi-squared = 22.006, df = 3, p-value = 6.505e-05

>

```

I then did Pairwise comparisons using Wilcoxon rank sum test to see which samples are different with each other.

```
> pairwise.wilcox.test(my_data$cell_type, my_data$group,
+                      p.adjust.method = "BH")

    Pairwise comparisons using Wilcoxon rank sum test

data:  my_data$cell_type and my_data$group

     CAF1   CAF2   NOF1
CAF2 0.0058 -      -
NOF1 0.2809 4e-05  -
NOF2 0.6818 0.0016 0.3811

P value adjustment method: BH
>

```

My Question is; how I can extract the genes behind the difference between each pair of samples? I mean something like differential expression by non parametric test. also I am not sure if I am right in using log cpm values for these test

If possible give me a hand to correct myself and get intuition

Thanks

# Answer
All of the genes are behind the difference, you're doing a rank test comparing all of them against each other at once. Your results have more to do with samples having slightly different outlier genes than any interesting biology. You don't have enough replicates to do a non-parametric test.
> 2 votes

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


title: Path to career in bioinformatics?

text:
# Question
Title: Path to career in bioinformatics?
I've been completing various courses on bioinformatics through edX and Coursera, and I've made a strong footing in the material.

I'm not sure what to pursue next in bioinformatics that would bring me closer to research in the field or employment.

I don't have bioinformatics at my university and thus am unsure of how to progress.

I have a background in mathematics, chemistry, physics and computer science. So, what would you recommend I do next?

# Answer
The Bioinformatics Stronghold at Rosalind is an excellent resource for building your understanding of many concepts and algorithms you'll see frequently in bioinformatics and genomics.

I'd also say that, if you can make it work, an internship where you can "get your hands dirty" with real problems is a very valuable experience. This option may not be available to everyone due to a variety of factors, but when it is an option it provides a tremendous opportunity to build experience with little risk to the intern or the employer.
> 3 votes

---
Tags: genomics
---


# Question
Title: Reads mapped to exonic, intronic and intergenic regions
After the alignment step I checked the rnaseq metrics of all the samples. Among 40 samples three samples show high percentage of reads mapped to intronic regions. What could be the reason?

```
Samples   Exonic            Intronic    Intergenic
Sample1 545479 (12.8%)  2512309 (58.8%) 1217201 (28.5%)
Sample2 372836 (8.3%)   2556934 (56.8%) 1573032 (34%)
Sample3 529618 (7.3%)   3934615 (54.5%) 2750720 (38.1%)

```

I also see that 35-40k genes among 58k genes having "zero" read counts. Is this due to contamination? What could be the reason? What does reads mapping to exonic, intronic, intergenic tell?

Update: I used hisat2 for alignment. I used human genome "grch38\_snp\_tran" from hisat2 website. Libraries were generated using ribosome depletion kit.

# Answer
> 3 votes
As @DevonRyan mentioned, it's very likely that those samples were degraded, which is good justification for excluding them from subsequent analysis.

# Answer
> 2 votes
I would expect at least 30% of reads from a total-cell, ribo-depleted RNA-seq to be exonic. Less suggests something when wrong.

As well as degradation, another explanation would be contamination with genomic DNA.

RNA obtained from the nuclear or cytoplasmic fraction of the cell might have a exonic content lower than 30%, while poly-A selected RNA would be expected to have higher, more like 66%.

---
Tags: rna-seq, sequence-alignment, rna, reads
---


# Question
Title: "Sequence Duplication Levels" module still fails after pre-processing Illumina data
I want to ask about why the sequence duplication levels are high after I trimmed by using Trimmomatic? I am using the following Trimmomatic operations: `HEADCROP = 19 TRAILING = 20 MINLEN = 66`.

How can i solve this problem? Thank You.

# Answer
To answer your direct question, there are a few reasons why there might be high levels of sequence duplication. From the FastQC help:

> *The underlying assumption of this module is of a diverse unenriched library. Any deviation from this assumption will naturally generate duplicates and can lead to warnings or errors from this module.*

* As @DevonRyan mentioned, with certain sequencing protocols such as RNA-Seq, two sequence reads at exactly the same location aren't that uncommon. This isn't a problem with RNA-Seq data, or with Trimmomatic, or with FastQC. It's just that this kind of data violates the assumption, and therefore should be ignored in those circumstances.
* PCR duplicates are another possible cause. PCR duplicates can give the false impression of high coverage at a particular locus when in fact it's just a single observed read that has been duplicated many times (see here for more details). PCR duplicates can usually be detected and removed if your analysis involves mapping to a reference genome. But whether this is actually a problem you need to fix depends on what type of data you have and what types of analysis you want to do.
* Large numbers of adapter dimers or rRNA may be present in your sample.

But I think it's also important to address how quality control (QC) is run. It can be tempting to run and re-run QC tools like Trimmomatic until all errors go away, but to be blunt these tools cannot think for you. For example, it's *possible* to get rid of most adapters by aggressively cropping/trimming both ends of each read, but you'll likely throw away a lot of good data that way. You may want to look into Trimmomatic's ILLUMINACLIP operation. It's also may be tempting to crop/trim reads aggressively if there are compositional biases near the beginning or end of the read. In fact, random hexamer priming can cause the **Per-base Sequence Content** module to fail on almost any RNA-Seq sample. Again, from the FastQC help (emphasis mine):

> *Libraries produced by priming using random hexamers (including nearly all RNA-Seq libraries) and those which were fragmented using transposases inherit an intrinsic bias in the positions at which reads start. This bias does not concern an absolute sequence, but instead provides enrichement of a number of different K-mers at the 5' end of the reads. Whilst this is a true technical bias, **it isn't something which can be corrected by trimming** and in most cases **doesn't seem to adversely affect the downstream analysis**. It will however produce a warning or error in this module.*

So in other words, it's best to determine what can cause each FastQC module to fail, investigate whether this is actually a problem for your data set (referring to documentation as needed), and make a deliberate QC plan that addresses the issues that need attention.
> 6 votes

# Answer
FastQC assumes that all samples are for whole genome sequencing and will flag them as failed if they differ too much from that assumption. This will, for example, cause essentially all RNA-seq, ChIP-seq, and ATAC-seq samples to fail in one module or another. This is not any cause for concern and is completely expected. Primarily concern yourself with whether all of your samples are similar in their metrics.
> 6 votes

---
Tags: illumina, data-preprocessing, trimming, fastqc
---


title: how to write italic script in Rstudio

text:
# Question
Title: how to write italic script in Rstudio
I have got the problem on how to write the italic script in RStudio. I used the script in the screenshot and got the error "Error in grobs\[\[i\]\]: subscript out of bounds."

# Answer
> 4 votes
The title is misleading as this error doesn't have anything to do with making the font italic. In Seurat v2, FeaturePlot does not return a ggplot2 object by default, so `p` in your case is `NULL`. You need to set `do.return` to `TRUE` in the FeaturePlot call.

You should instead do:

```
library(Seurat)

p <- FeaturePlot(pbmc_small, head(pbmc_small@var.genes), do.return = TRUE)

for(i in 1:length(p)) {
  p[[i]] <- p[[i]] + theme(plot.title = element_text(face = 'italic'))
}

cowplot::plot_grid(plotlist = p)

```

---
Tags: rna-seq, single-cell
---


# Question
Title: Drug with as many targets as possible through a system biology graph?
I have 4 pandas data frames containing the following data:

1. affinity\_score1: affinity between each couple (Drug-Target) measured using tool1
2. affinity\_score2: affinity between each couple (Drug-Target) measured using tool2

The two tools (tool1, tool2) give measure in different scales and are not 100% accurate.

3. Drug\_similarity: a similarity matrix between drugs
4. Target\_similarity: a similarity matrix between targets.

Here also the tool used to measure the similarity is not 100% accurate.

```
import pandas as pd

affinity_score1 = [{'Drug0': -24.5, 'Drug1': -17.1, 'Drug2': -13.1},
           {'Drug0': -20.0, 'Drug1': -21.6, 'Drug2': -18.1},
           {'Drug0': -19.0, 'Drug1': -15.4, 'Drug2': -20.0}]

affinity_score2 = [{'Drug0': 1.5, 'Drug1': 0.8, 'Drug2': 0.4},
           {'Drug0': 1.2, 'Drug1': 1.3, 'Drug2': 0.9},
           {'Drug0': 0.8, 'Drug1': 0.7, 'Drug2': 0.9}]
affinity_score1_df = pd.DataFrame(affinity_score1).rename(index={0:'TargetX',1:'TargetY', 2:'TargetZ'})
affinity_score2_df = pd.DataFrame(affinity_score2).rename(index={0:'TargetX',1:'TargetY', 2:'TargetZ'})
Drug_similarity = [{'Drug0': 1, 'Drug1': 0.8, 'Drug2': 0.2},
           {'Drug0': 0.8, 'Drug1': 1, 'Drug2': 0.9},
           {'Drug0': 0.2, 'Drug1': 0.9, 'Drug2': 1}]
Target_similarity = [{'TargetX': 1, 'TargetY': 0.8, 'TargetZ': 0.2},
           {'TargetX': 0.8, 'TargetY': 1, 'TargetZ': 0.9},
           {'TargetX': 0.2, 'TargetY': 0.9, 'TargetZ': 1}]
Drug_similarity = pd.DataFrame(Drug_similarity).rename(index={0:'Drug0',1:'Drug1', 2:'Drug2'})
Target_similarity = pd.DataFrame(Target_similarity).rename(index={0:'TargetX',1:'TargetY', 2:'TargetZ'})

```

**I am trying to find the Drug who can have as many Targets as possible.**

So far, a simple approach could be to select the Drug having the highest average score in both affinity scoring tools.

This approach does not use the information in the similarity matrices. Given that if a Drug1 is successful with TargetX, and TargetX is similar to TargetY then Drug1 can also easily be successful with TargetY. Considering the inaccuracy of the scoring tools, I think including this information could improve the accuracy of the answer, eliminating some false positives. Also, if a Drug from a cluster of similar drugs (information that can be derived from the drugs similarity matrix) is the only one predicted to have good affinities with some targets, then this is more likely to be a false positive. But I don't know how to efficiently integrate all the information. Scipy may help with clustering the two similarity matrices.

Another way of analysis is to set affinity and similarity thresholds, in which case each data frame could be transformed into a network:

* Each node is a drug or a target.
* Two types of edges: similarity (between drugs or targets) and affinity (drug-target)
* The connected nodes have passed the thresholds.
* The edges have weights: affinity or similarity score.

The problem then becomes finding the most connected node (a drug node) with the highest weight values (scores on tool1 and tool2). Connected nodes through similarity may then be used to validate the number of targets for each drug. It seems a more complicated approach but more rigorous with most likely higher positive predictive value.

Any python package to implement such an approach? And also a visualization tool.

The original data has 38 Targets with 806 drugs.

# Answer
I don't know any python (or R package) to do this), but this is how I would do it based on these matrices:

1. Do a multidimensional scaling (MDS) plot of the target\_similarity matrix. Explore the visualization:
2. Find which drug targets the drug closest to the 0,0 point in the 2D space

However this doesn't say anything about the affinity between the drug and the target. If you wanted to add this then you could weight how much close is the distance between the drug and the target vs how close are they to the center of the MDS (a simple division could work if the distance is scaled between 0 and 1 and averaged).

Note that depending on the original information used to calculate the similarities and the affinities it could be done in other ways
> 2 votes

---
Tags: python, drugs
---


# Question
Title: How to plot genomic.fna fasta file in R using Gviz?
I downloaded the genome of Staph aureus as DNA sequence from NCBI. I would like to visualize it using Gviz. Is Gviz the right tool to visualize genomes? What are the necessary steps from the sequence file to an object that can be used in Gviz? I searched the user guide and the reference manual. It is stated that it is possible to use fasta file, but not how to use them.

We will generate lots of genomes that we will have to annotate and compare. S aureus is not S aureus. There is a spectrum of phenotypes, they all have a different genomes. And we will also not be limited to SA. The question here is not how to annotate the genome, but rather the steps necessary to visualize something and get going. It is my first steps with genomic data, I don't even know what makes sense.

Here is the head of the file:

```
>NC_007795.1 Staphylococcus aureus subsp. aureus NCTC 8325 chromosome, complete genome
CGATTAAAGATAGAAATACACGATGCGAGCAATCAAATTTCATAACATCACCATGAGTTTGGTCCGAAGCATGAGTGTTT
ACAATGTTCGAACACCTTATACAGTTCTTATACATACTTTATAAATTATTTCCCAAACTGTTTTGATACACTCACTAACA
GATACTCTATAGAAGGAAAAGTTATCCACTTATGCACATTTATAGTTTTCAGAATTGTGGATAATTAGAAATTACACACA
AAGTTATACTATTTTTAGCAACATATTCACAGGTATTTGACATATAGAGAACTGAAAAAGTATAATTGTGTGGATAAGTC
GTCCAACTCATGATTTTATAAGGATTTATTTATTGATTTTTACATAAAAATACTGTGCATAACTAATAAGCAAGATAAAG
TTATCCACCGATTGTTATTAACTTGTGGATAATTATTAACATGGTGTGTTTAGAAGTTATCCACGGCTGTTATTTTTGTG
TATAACTTAAAAATTTAAGAAAGATGGAGTAAATTTATGTCGGAAAAAGAAATTTGGGAAAAAGTGCTTGAAATTGCTCA
AGAAAAATTATCAGCTGTAAGTTACTCAACTTTCCTAAAAGATACTGAGCTTTACACGATTAAAGATGGTGAAGCTATCG
TATTATCGAGTATTCCTTTTAATGCAAATTGGTTAAATCAACAATATGCTGAAATTATCCAAGCAATCTTATTTGATGTT

```

# Answer
I am not familiar with Gviz, but I am familiar with a variety of tools for visualizing genome features. At the scale of entire microbial genomes, the primary sequence itself is an overwhelming amount of information and can't really be visualized in a meaningful way by itself. I mean, you could color the nucleotides, but few patterns in DNA can be discerned by an unassisted human eye.

What we're typically interested in are features encoded in the genome sequence: protein-coding genes, non-coding RNAs, transposable elements, regulatory elements, binding sites, epigenetic modifications, and so on. A particular scientist is usually only interested in a subset of these features, based on what particular questions they are investigating. In your comment you say you "*would like to annotate protein coding regions and compare them to other genomes*", so it seems like we're on the same page here.

Visualization tools often display the genome linearly (even, for example, if it's a circular bacterial chromosome), zoomed out so that the nucleotide sequence isn't visible, and overlayed with boxes, arrows, or other glyphs representing the genes or other features of interest. Here are some examples:

In each of these cases, the visualization software is not doing the work to determine the positions of the genes and other genome features. This was done independently using tools designed for this task, storing the information in an **annotation** file in GFF3, GTF, or BED format. The visualization tools need both the genome sequence (in Fasta format) and one or more annotation files (in GFF3, GTF, or BED format) to make meaningful visualizations.
> 2 votes

---
Tags: r, genome, visualization, gviz
---


title: Genomics Statistics Problem

text:
# Question
Title: Genomics Statistics Problem
I have used a python script to identify target sequences in a DNA sequence file. There are two classes of sequence: coding and non-coding.

I have identified 728 sequences of interest. 597 of these fall into the coding regions and 131 of these fall into the non-coding regions. This is the equivalent of 18% non-coding. The total non-coding region in the sequence file is 13%.

Is there a statistical tool to demonstrate the python script identified target sequences in a non-random fashion?
If the script identified sequences that were randomly distributed then 13% of them would have been found in the non-coding region, from a total of 728 this seems like it should be reliable.

# Answer
> 3 votes
If your null hypothesis is that the event of being selected by your program is statistically independent from the event of being non-coding, you can use Fisher's exact test to reject the null hypothesis.

Suppose you have 10,000 (I'm making that up but presumably you have the real number) total sequences, 8,700 coding and 1,300 noncoding. Given that your program selects 728 sequences, you would expect around 13%\*728 ~= 95 noncoding sequences. But your program picked 131 noncoding sequences and you want to know if 131 is far enough away from 95 to conclude that the difference is not just due to chance.

The probability of selecting 131 or more of these noncoding sequences (given that we are choosing 728 total sequences) is given by the so-called hypergeometric distribution. The two-sided test will also incorporate the probably of getting an equally unlikely low number (in this case, somewhere around 57) or less.

In R for example, you can calculate this using the 2-by-2 contingency table as follows:

```
fisher.test(rbind(c(597, 131), c(8103, 1169)))

```

which gives:

```
p-value = 7.399e-05

```

So in this situation (and assuming the probability model described), there's around a 0.0074% chance of observing a ratio as "extreme" as the one you observed, which indicates it's unlikely that the two events are independent.

Note: As others have mentioned, a binomial test could be appropriate too. If your number of total sequences is very large, then the two results will converge. But if it isn't, then Fisher's exact test will be more appropriate (essentially since you're sampling without replacement).

# Answer
> 2 votes
A Poisson or binomial would suffice.

I'm not at all sure the difference will significantly deviate from randomness.

Perl: 95% confidence interval of a non-parametric bootstrap for e.g. 1000 replicates is an alternative approach. I'd be more comfortable coding that solution in Perl than I would in Python, so is unlikely to be of much use here. Anyway, I would use a 2-dimensional array for your total data-set, the second dimension being coding/non-coding, I would then use the inbuilt random-number generator and multiply by the size of the array to randomly pluck out each element. 95% CI of the mean works by ranking the data and removing the top and bottom 2.5% and taking the interval for a single replicate (which is repeated 999 times more).

There will be a deep learning approach to this data, which would be easier to implement in Python, but I'd need to think about this because it is not simple classification.

Again my warning is that given the sample size, the result at best might be borderline. There are likely to be lots of interesting in your data, but coding/non-coding may not be one of them.

---
Tags: phylogenetics, phylogeny, statistics, genomics
---


title: Changing this function to work

text:
# Question
Title: Changing this function to work
I have this matrix of expression data

```
> head(data[,1:6])
  Symbol      cl1    cl1.1     cl1.2     cl1.3    cl1.4
1   ACTB 14.07922 14.25211 13.547206 14.046814 15.33551
2 ATP5F1 10.43386 10.34219  9.671036  9.862641 10.39786
3   DDX5 13.99250 13.29734 13.706894 13.084000 13.73379
4  EEF1G 14.35303 14.26112 13.523434 13.327333 13.97747
5  GAPDH 16.07695 15.61435 15.454911 15.233670 15.90100
6    NCL 14.35149 13.32919 13.051861 12.704973 13.55680
>

```

I have 24 samples as cl1 and 30 samples as cl2

I want to do a t.tset, I have this function but I am not sure Why I am getting error by

```
volcan <- function(data, cl1, cl2) {
   tt <- t(apply(data[, -1], 1, function(d){z <- t.test(d[cl1], d[cl2]); c(z$estimate, z$p.value)}))
     tt[, 1] <- tt[, 1] - tt[, 2]
     tt <- tt[, -2]
     tt <- cbind(symbol=data[, 1], as.data.frame(tt), rank=NA)
     tt[, 3] <- -log(tt[, 3], 10)
     colnames(tt)[2:3] <- c("log.ratio", "neg.log(t.test)")
     tt[, 4] <- sqrt(tt[, 2]^4  tt[, 3]^2)
     tt <- tt[order(tt[, 4], decreasing=T), ]
     tt
}
    Error: unexpected symbol in "tt[, 4] <- sqrt(tt[, 2]^4 tt"

```

When I am running function as whole says that

```
Error: unexpected '}' in "}"

```

When I am putting

> tt= volcan(data,cl1,cl2)

and doing line by line going down after

```
tt[, 4] <- sqrt(tt[, 2]^4  tt[, 3]^2)

```

Says

```
Error: unexpected symbol in "tt[, 4] <- sqrt(tt[, 2]^4  tt"

```

I am not sure why, I put comma in

```
 tt[, 4] <- sqrt(tt[, 2]^4 , tt[, 3]^2)

```

But another error

# Answer
> 5 votes
EDIT: I'm editing my answer as I don't have comment privileges and don't wish to add multiple answers.

Please don't change the entire premise of the question after you get an answer. If you do, retain the original content so the post makes *some* sense to people.

Run `?sqrt` and check how the function is to be invoked, and how you're invoking it. Does it even accept 2 input parameters?

---

## Original Content

You are subsetting the 2D object `data` improperly. How can data be subset as both `data[,1:6]` and `d[data[,1:24]]` where `d` = `data`? Fix that and your r code will work OK (syntactically at least).

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


# Question
Title: Visualising gene expression across cell type and conditions in one plot, in Single Cell Sequencing data
I am new to single cell sequencing data analysis but I have basic programming background in R and python. I want to be able to plot differential expression of two genes, Gene1 and Gene2, across three cell types and across three conditions on one plot. I have followed this tutorial to perform QC and various differential analyses using Seurat package on my data up till now.

My data set is one object, containing information about the three cell types under three various conditions. The metadata slot of my data set contains information about my cell types as well as the conditions under which they are tested. Using the following DotPlot commands I am able to generate separate plots of gene expression with respect to cell type and with respect to condition:

```
Seurat::DotPlot(object = object, genes.plot = c("Gene1", "Gene2"), plot.legend = TRUE, group.by = "cell_types")

Seurat::DotPlot(object = object,
        genes.plot = c("Gene1", "Gene2"), plot.legend = TRUE, group.by = "tissue_types" )

```

For elegant representation of my data, I want to be able to have the information of the following plots into one plot.

I will be very grateful for any hint on how I can proceed to achieve an elegant presentation as well as information about other possible ways to graphically represent gene expression data for publication or presentation purposes.

Thank you in advance for your help.

# Answer
> 2 votes
Try using `SplitDotPlotGG`, this will allow you to color the dots by some additional variable (eg condition).

In your case you would want to do something like this:

```
library(Seurat)

SplitDotPlotGG(object = object,
               genes.plot = c("Gene1", "Gene2"),
               plot.legend = TRUE,
               group.by = "cell_types",
               grouping.var = "tissue_types")

```

This would give a result similar to below (an example from https://satijalab.org/seurat/immune\_alignment.html).

---
Tags: scrnaseq, seurat
---


title: SNP vs Point Mutation

text:
# Question
Title: SNP vs Point Mutation
What is the difference between a Single Nucleotide Polymorphism (SNP) and a point mutation? I am quite confused in understanding these term as both of them refer to one base difference from the reference sequence?

# Answer
The difference between the two depends on to whom you talk ;)

You are right: both refer to one base difference from the sequence. Usually the term "mutation" is used if the change has an impact on the phenotype. The "P" in "SNP" means "polymorphism". A lot of people use this term if more than 1% of the individual in a given population have this variant.

My recommendation: only use "mutation" if you are sure that your variant has an important impact on the phenotype. Only use "SNP" if you are sure that it is a commonly occurring variant. In any other case just use "variant". It is slowly becoming more common to talk about "SNV" (Single Nucleotide Variant) than about "SNP".
> 6 votes

---
Tags: snp, terminology, variation
---


# Question
Title: log rank p value interpretation from gepia 2
How do i interpret the overall survival when the p value reported is not significant going by the standard convention .

Not sure what kind of test it runs while comparing the expression may be one non-parametric test , but the survival is kaplan meier on cox ph model which they have mentioned .

So is it a significant one if i have to tell in context of survival ?

# Answer
> 1 votes
The plot on the left shows that the gene has very different expression between two groups. However, according to the graph on the right, these high expression is not directly related to survival because both groups have similar survival time.
The gene expression so high might be to compensate for other deficiencies or problems.

---
Tags: rna-seq, figure-explanation
---


# Question
Title: What is the difference between NGS, 2GS, SBS and HTS?
I've come across a bit of confusion about the initialism NGS, so think it would be a good idea to clarify this term (and similar terms like 2GS, SBS, and HTS) for this site with a bit of discussion.

What are the most common initialisms used to describe different sequencing technologies?

# Answer
> 5 votes
Here are my attempts at definitions:

Sanger: A method of sequencing that depends on chain-terminatiing dideoxynucleotides. This sequencing uses the differential flow of DNA sequences of different lengths through a gel to determine the original DNA sequence, producing a single sequence per reaction container.

NGS: Next-generation sequencing, also referred to as 2GS (second-generation sequencing). This term is used to describe the first wave of sequencing technologies that followed Sanger sequencing technology. The use of NGS has become more confusing with the advent of long-read sequencing, because it's a common assumption that "next-generation" refers to the most recent technology (which is incorrect in this case).

HTS: High-throughput sequencing. This term describes any type of sequencing technology that produces large amounts of data, usually in the form of millions of different sequences produced from the same sequencing run.

SBS: Sequencing by synthesis. This term describes a method of sequencing that depends on the synthesis of \[DNA\] bases in order for sequencing to be carried out. This definition can extend into long-read sequencing (e.g. PacBio sequencers depend on synthesis during sequencing), but is more typically associated with only the second-generation sequencing technology.

# Answer
> 5 votes
I am not sure how appropriate it is at this point to still refer to sequencing as next-generation sequencing. The leading NGS technology is Illumina/Solexa that has been around for over 10 years at this point. 454 was around even earlier. It's not really "next" at this point.

Opinions aside, I would also add "third-generation sequencing" to that list, referring to long read technologies like PacBio and Oxford Nanopore. See this question on Biology SE for more details.

# Answer
> 0 votes
The terms "second-generation sequencing", "third-generation sequencing", "short-read sequencing", and "long-read sequencing" are used commonly but (as far as I know) lack commonly used acronyms/initialisms.

In the forensic genetics community, the acronym **MPS** is used commonly to refer to "massively parallel sequencing." This definitely includes 2G/short-read sequencing (such as Illumina). It's less clear whether MPS includes 3G/long-read sequencing such as PacBio or Oxford Nanopore, but since these are rarely used in a forensics context it's a moot point.

Of course, it's important to acknowledge (as @burger has done) that there is a lot of imprecision and opinion regarding the appropriate scope of each term. Some folks seem to want a term to refer to everything post-Sanger, which includes an incredible diversity of platforms with respect to sequencing technology, read length, throughput, and error profiles. Some very reasonable questions should make us think seriously about our chosen terminology: How long do (for example) Illumina reads have to get before they're no longer classified as short reads? How much throughput is required of a sequencer to be classified as "high" throughput? When do these distinctions *really* matter?

---
Tags: ngs, terminology, hts
---


# Question
Title: Remapping genomic coordinates to account for indels
I'm interested in obtaining coding sequences of my favourite gene in all individuals from the 1000Genomes (and similar projects). I use GATK to get the right subset of variants, vcf-consensus to map these variants onto the reference genome and finally samtools to extract the individual exons. This works fine if the variants are SNPs but if there are any indels, this changes the coordinates of exons and I end up getting the wrong region. Is there any generic way of remapping genomic coordinates to account for the changes created by indels?

# Answer
> 8 votes
I think that you need a LiftOver Chain file to transform your coordinates. You can obtain such a file using `bcftools consensus` with the `-c` parameter:

```
-c, --chain <file>         write a chain file for liftover

```

Then you can use it to transform coordinates in various genomic formats using CrossMap.

# Answer
> 2 votes
If the goal is to extract the consensus sequence for given regions, this is how it works in 2019.

`bgzip` and index your `vcf` file.

```
$ bgzip -c input.vcf > input.vcf.gz
$ tabix input.vcf.gz

```

Create a `regions.txt` which contains one region per line in the format `chr:from-to`. If you have already a `bed` file, you can use this little `awk`script to create it:

```
$ awk '{print $1":"$2+1"-"$3}' input.bed > regions.txt

```

Now run this combination of `samtools` and `bcftools`:

```
$ samtools faidx -r regions.txt genome.fa|bcftools consensus input.vcf.gz -o consensus.fa

```

---
Tags: variation, structural-variation
---


# Question
Title: How to manipulate a reference FASTA or bam to include variants from a VCF?
I have some software which takes fastas as the input. I need to include SNVs and InDels from a VCF into the reference hg38 and then use this.

The problem is, I don't know of an algorithmically sound way to do this.

1. Are there any existing software packages which could do this efficiently? Is it easier to output a FASTA, or a bam (and then convert to a FASTA)?
2. What about if I wanted to do the same with a bedpe of germline structural variants?

# Answer
You could convert VCF to BED via `vcf2bed --snvs`, `vcf2bed --insertions`, and `vcf2bed --deletions`, and then use `samtools faidx` by way of a wrapper script to convert BED to FASTA, e.g.:

```
$ vcf2bed --snvs < variants.vcf | bed2faidxsta.pl > snvs.fa
$ vcf2bed --insertions < variants.vcf | bed2faidxsta.pl > insertions.fa
$ vcf2bed --deletions < variants.vcf | bed2faidxsta.pl > deletions.fa

```

You need FASTA files for your reference genome, which have been indexed with `samtools faidx`, e.g., for `hg38`:

```
$ cd /foo/bar/baz
$ wget ftp://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/*.fa.gz
$ for fn in `ls *.fa.gz`; do gunzip $fn; done
$ for fn in `ls *.fa`; do samtools faidx $fn; done

```

Once you have indexed FASTA files somewhere on your file system, you can pipe BED to the `bed2faidxsta.pl` script, to get out FASTA sequences.
> 4 votes

# Answer
GATK has a solution that might work for you: FastaAlternateReferenceMaker, which : "Given a variant callset, this tool replaces the reference bases at variation sites with the bases supplied in the corresponding callset records."

> **Input**
>
> The reference, requested intervals, and any number of variant ROD files.
>
> **Output**
>
> A FASTA file representing the requested intervals.
>
> **Usage example**
>
> ```
> java -jar GenomeAnalysisTK.jar \
>    -T FastaAlternateReferenceMaker \
>    -R reference.fasta \
>    -o output.fasta \
>    -L input.intervals \
>    -V input.vcf \
>
> ```
> 7 votes

# Answer
There's a `vcf2fq` sub-program that was written as part of vcfutils to convert a VCF file into a fastq file given a reference sequence. Unfortunately this doesn't work properly with INDELs (it will just mask them, rather than actually converting them), so I wrote a modification to implement INDEL correction as well:

```
./vcf2fq.pl -f <input.fasta> <all-site.vcf> > <output.fastq>

```

It adds and removes sequence from the reference, as specified by the INDEL information in the VCF file. The code keeps a record of where in the reference sequence the INDELs should go, then inserts (or deletes) the INDELs, updating the reference position as necessary. It should work with things like deletions right next to insertions, but hasn't been extensively tested.
> 3 votes

# Answer
Here's an update how this work in 2019. (The answer is quite similar to my update in this thread):

`bgzip` and index your `vcf` file.

```
$ bgzip -c input.vcf > input.vcf.gz
$ tabix input.vcf.gz

```

Now use bcftools consensus:

```
$ bcftools consensus -f genome.fa input.vcf.gz > consensus.fa

```
> 2 votes

---
Tags: sam, fasta, vcf, reference-genome, variation
---


# Question
Title: Removing genes with less than a correlation cut-off between two matrices
I have two matrices like this:

```
> head(DEG_log[,1:4])
            cl1       cl1       cl1       cl1
ALB    4.653796  7.046317  7.732135  7.758708
AQP9   7.798208  8.231112  8.734799  8.054087
CALML5 5.791524  9.058341  6.977369  9.503477
CCL4   9.739961  7.663303  8.934939  7.324958
CDK6   9.832115  9.710043 11.439345 10.330253
CHGA   5.139340 12.024606  5.988773  8.696836
>

```

I have the same patients and genes in two separate matrices coming from two sequencing panels with the same chemistry. Each sample has a match in another panel. I want to remove genes that have low correlation between two matched samples.

I have normalised and log-transformed the data by DESeq2.

I just plotted the correlation between each pairs of samples:

```
library("ggpubr")
ggscatter(my_data, x = "sample1", y = "sample1_1",
          add = "reg.line", conf.int = TRUE,
          cor.coef = TRUE, cor.method = "pearson",
          xlab = "onco", ylab = "biomarker)")

```

Like this:

by:

```
corr_test(a, b, 0.95)

```

I can calculate correlation, but I don't know how to remove genes with low or non-significant correlation between matched samples.

# Answer
Nice data and good question.

It is called "outlier" analysis or the analysis of the residual. Its an easy analysis, but . You calculated the standard deviation using the regression slope as the 'mean' and retain data within 1.96 standard deviations. You can plot the data like this and if you are lucky will see the residual forming a nice 'normal distribution'. Outliers outside 1.96 sd can be formally excluded.

Calling (and plotting) the residual from the regression slope will be possible in R, but I use external packages because its easy. BTW you should also try Spearman's coefficient to check you are getting the same answer.
> 3 votes

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


# Question
Title: How is consensus alignment for OLC assembly usually implemented?
## Background

In C++ I want to implement a simple class that performs overlap layout consensus assembly but I can not figure out the most logical data structure to use for the consensus alignment step.

This is the step where all of the reads from the overlap graph are aligned and the consensus is called (slide 48). An example alignment is below. Asterisks above the column indicate some disagreement among all reads. Alignment position is shown above and below for convenience.

```
                         *    * *         *
     pos           111111111122222222223333333333444
         0123456789012345678901234567890123456789012
read 1   TAGATTACACAGATTACTGAT TGA
read 2      ATTACACAGATTACTGACTTGATG CGTAAACTA
read 3                   tTGACTTcATGGCGTAA CTA
read 4                     GACTTGATGGCGTAA CTAAG
read 5                           ATGGCGTAA CTAAGATCG
     pos           111111111122222222223333333333444
         0123456789012345678901234567890123456789012

```

## Question:

**Q1:** What is the most logical way to implement the consensus alignment step in a OLC assembler?

**Q2:** What are the pros and cons of this and other possibilities?

## Ideas

### Idea - Vector of vectors. One column per position.

My initial idea for a data structure to store the alignment is to use one vector per read to store the read's character at that position. We would end up with j arrays (number of reads) all of length i (the total length of the overlapped alignments). This, however, would be space-inefficient if the locus were large and many positions in the vectors would contain a null value. In the example alignment above the first 23 positions of read 5 would be null values.

### Idea - Graph. One node per position.

In another implementation idea, each position in the alignment would be represented by a node in a directed graph. Each node would have a property keeping track of the nucleotide counts at each position and might keep the read id's stored as well if there is a good reason to do so.

# Answer
> 2 votes
For short reads, mainstream assemblers don't have a consensus step. The contig path directly spells the consensus sequence.

For long reads, the most popular way to derive a consensus is via partial-order alignment (POA). This is used by error correctors (similar to consensus) like pbdagcon and falcon\_sense, and consensus tools like canu, racon, wtdbg2 and nanopolish. Basically, you start with an empty directed acyclic graph (DAG), align each read in turn to the graph and augment the DAG by inserting unaligned segments. When you align all reads onto the DAG, the most weighted path is taken as the consensus. You should read Jared Simpson's blog on this topic. It is much more detailed than "basically". You should also read the papers on all the mentioned tools above if you want to implement by yourself.

An alternative approach is Gene Myers' realigner (PDF here). It is not often used, though.

---
Tags: assembly, algorithms, c++
---


# Question
Title: Calculating Z-score from logCPM values using edgeR
I have the raw counts for RNA-Seq data. I converted counts data to logCPM using edgeR package.

Lets say I have a dataframe `A` with 15000 genes as rows and 100 samples as columns with counts data. I have used following code from edgeR tutorial.

```
logCPM <- cpm(A, prior.count=2, log=TRUE)
head(logCPM)[1:5,1:5]
               Sample1          Sample2          Sample3         Sample4
A1BG          -1.019429        0.8333739        -1.192960       -0.6717287
A2M            7.622280        7.8057767         7.932570        8.6409077
A2ML1          1.004541        1.0909119         5.462613        6.4149348
A4GALT         2.126502        1.8403954         5.322601        2.9469345
AAAS           5.785888        5.1639501         5.259789        5.3248026
               Sample5
A1BG          0.3121975
A2M           9.1949911
A2ML1         9.4740594
A4GALT        2.9214979
AAAS          4.3708063

```

After converting counts to `logCPM` I see there is another command for scaling which is used for heatmaps.

```
logCPM <- t(scale(t(logCPM)))

```

It is given there that `scaling each row (each gene) to have mean zero and standard deviation one` in page 15 edgeR

Is this nothing but z-score? or Is there any otherway to calculate z-score from logCPM data?

thanq

# Answer
> 1 votes
Yes this is nothing but a z-score. The scale function calculates column z-scores, but since you use t (transpose), it is actually performed on rows. There are many other ways to calculate z-scores, but if done correctly they will all give the same results.

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


title: Non-random access on a fastq file

text:
# Question
Title: Non-random access on a fastq file
There was a nice question on random access on a fastq file, however I have the opposite problem. I want to programatically access non-random fastq entries quickly (using C++).

Say I have several pairs of `fastq.gz` files that could be hundreds of millions or \> 1 billion read pairs total. I know the number of read pairs in each file and I have a list of 10-20 reads that I want to access from each file. I only know the 0-based index of the fastq entry that I want to access, not the read ID. So if I want the 234th and 863rd reads the naive solution is to open the fastq file and scan through 864 records. I need to perform some tasks with these reads, then I move on and access/process another batch of reads.

* I need to perform the above task tens of thousands of times, so the options are:
    + 1. **Low RAM usage. Many file IO ops. No disk usage:** Opening/scanning/closing the `fastq.gz` files then processing the desired reads **tens of thousands of times.**
    + 2. **Extreme RAM usage. One round of file reading. No disk usage:** Read in each `fastq.gz` file once, store all of the desired reads in RAM and access and process those reads as needed
    + 3. **Low RAM usage. One round of file reading. Some disk usage:** Coming up with an intermediate file format similar to a `bam` that will allow me to parse and sort the reads with their groups and other metadata later on.

If it makes a difference the fastq.gz files are on spinning disks, not SSDs.

# Answer
> 2 votes
You don't have a RAM problem because you only have about:

> 10-20 reads that I want to access from each file.

Please go for (2). If you are concerned about RAM usage, you may selectively write some of your data back to disk.

Tens of thousands of 10-20 reads is not a problem for modern computing.

---
Tags: fastq, algorithms, c++
---


# Question
Title: biological fluorescence supervised machine learing
Some general questions...

So it appears with most classification learners one must come up with a series of quantifiable variables to associate with each observation. This might include mean and std deviation intensity values, texture, and zernike polynomials. Does anyone have any other suggestions as to how one might quantify 3-dimensional fluorescence pixel data features?

Once appropriate quantification's have been acquired for a given class one can then train a classifier. My next question is - does one then have to segment test images and create a matrix of the same feature variables from that segmented image to use the classifier? Or can one use the raw image itself?

What supervised machine learning techniques will evaluate the raw image data itself - either for training or testing? Can manually annotated data sets or clustered data sets provide a direct model for training data? That is to say, are there supervised machine learning algorithms which accept pure pixel data as models for training?

Are there reliable pre-trained classifiers that are already implemented for these type of fluorescence data sets?

Any resources or advice is greatly appreciated.

Thanks in advance for any input!

# Answer
> 2 votes
Whether to quantify things at all with depend on the technique you want to use. An end-to-end CNN (convolutional neural network) for example would just be fed the 3D image for training and predictions. However the major downside to CNNs is that you need a LOT of data to train them, though perhaps you can get lucky and find an existing CNN that you can use and just retrain the last couple layers.

Other than that, you're going to need to come up with as many (hopefully useful) metrics as you can. There are then a lot of techniques one could try, ranging from simple things like regression to random forests to SVMs, among many others. In general I think simple regression and random forests are the easiest to get started with.

# Answer
> 2 votes
Adding to Devon's correct answer: Depending on the amount of images you got, there is also the possibility to avoid the feature engineering by building a two level architecture:

* The first level identifies the regions of interest in your image and scales these regions to the same size.
* In the second level, the recognized areas are fed to a (multilayer) neural network/deep neural network to be classified.

By using this approach, the number of needed datapoints is lower compared to a plain CNN method.

The drawback is that identifying the regions of interest itself might be a difficult task as you have either to label these regions manually to apply a machine learning approach at this first level, or you have to "engineer" criteria which make the extraction possible.

I have no idea how these images look like, but maybe you are able to extract the regions using OpenCV.

# Answer
> 0 votes
> What supervised machine learning techniques will evaluate the raw image data itself - either for training or testing?

1. You can use convolutional neural networks (CNN). They will directly use the raw image data itself.
2. After applying CNNs, you can also use recurrent neural layers (RNN), to achieve better results (e.g. gated recurrent layers).

Hope this helps. As a reference you have this article:

> Gudla PR1,2, Nakayama K2,3, Pegoraro G1,2, Misteli T2. Article: SpotLearn: Convolutional Neural Network for Detection of Fluorescence In Situ Hybridization (FISH) Signals in High-Throughput Imaging Approaches. Send to Cold Spring Harb Symp Quant Biol. 2017;82:57-70. doi: 10.1101/sqb.2017.82.033761. Epub 2017 Nov 28.

---
Tags: 3d-structure
---


# Question
Title: How to set the position of groups in a Seurat object on a FeatureHeatmap plot
I am analysing singe cell sequence data and I have followed this tutorial, https://satijalab.org/seurat/pbmc3k\_tutorial.html to perform QC and various differential analyses using the Seurat package on my data up till now.

How can I rearrange on the resulting plot, the order in which the groups of my samples appear when I use the `FeatureHeatmap()` function to visualise the differential gene expression of Gene1 and Gene2 across the cell clusters?

```
FeatureHeatmap(object = object, features.plot = c("Gene1", "Gene2"),
               group.by = "Conditions",
               pt.size = 0.25, key.position = "top",
               max.exp = 3)

```

produces:

but I want:

# Answer
I don't think this is possible in Seurat v2, but in v3 you can change the factor levels of the grouping variable to change the plot order:

```
library(Seurat)

FeaturePlot(object = pbmc_small,
            features = head(VariableFeatures(pbmc_small), 2),
            split.by = 'groups')

```

Change the order:

```
pbmc_small$groups <- factor(pbmc_small$groups, levels = c('g2', 'g1'))

FeaturePlot(object = pbmc_small,
            features = head(VariableFeatures(pbmc_small), 2),
            split.by = 'groups')

```
> 5 votes

---
Tags: scrnaseq, seurat, ggplot2
---


# Question
Title: How do I carry out an ancestry/admixture test on a single VCF file?
*This is a question from /u/beneficii9 on reddit. The original post can be found here.*

> Through the Personal Genome Project, I have had my whole genome sequenced by Veritas, and have it in the form of a single VCF file for the whole genome and one BAS file for each chromosome. The reference genome associated with the VCF file is hg19. It has been helpful in health data; for example, I discovered I'm homozygous for the non-functional variant CYP-2D6 gene (rs3892097), which can render several common medications useless, and helps explain why some medicines didn't really work for me. My doctor has found this information very helpful.
>
> Unfortunately, I can't find any way of looking at admixture or ancestry. I've tried setting everything up using a combination of VCFTools, Plink1.9, and ADMIXTURE, but I can't get it to work. I think for ADMIXTURE you have to have a bunch of genomes sorted by geographical origin to compare your genome against, but I'm not sure how to do that, and what's online isn't very clear to me. So scratch that one off.
>
> I've tried converting the file to 23andme format (and at this /u/psychosomaticism has been very helpful). I did that (though it seems there were problems because of the way the VCF file was set up). But the websites that take the data want you to point them to your 23andme account, and that doesn't really work if you only have the file. 23andme doesn't provide for people who had their whole genomes sequenced. They want you to give them a saliva sample like everyone else.
>
> So, what can I do?

# Answer
> 5 votes
A modified implementation of Vivek's answer.

peddy is a Python package that samples an input `.vcf` at ~25000 sites and projects onto a principal component space built on 2504 thousand genome samples. The author has extensive documentation of the tool's features and a link to the preprint.

I downloaded the `.vcf` and `.vcf.tbi` for the NA12878 sample from Genome in a Bottle's ftp here. Then, created a custom `.ped` file, `NA12878.ped`, with the contents:

`NA12878 HG001 0 0 2 0`

At the command line:

`$ peddy --plot --prefix myvcf HG001_GRCh37_GIAB_highconf_CG-IllFB-IllGATKHC-Ion-10X-SOLID_CHROM1-X_v.3.3.2_highconf_PGandRTGphasetransfer.vcf.gz NA12878.ped`

The output files all have the prefix `myvcf.`, here's `myvcf.pca_check.png`

# Answer
> 4 votes
A suggestion from the UK Biobank QC methods white paper:

1. Create a set of SNPs common to your VCF and the 1000 genomes phase 3 call set.
2. Perform PCA of the 1000 genomes samples using eigenstrat smartpca. You might have to convert to binary plink format.
3. Project your genotypes onto this pre-computed PCA space and visualize using ggplot to see which cluster you fall in.

# Answer
> 2 votes
Ancestry testing is a tricky subject. I spent a good chunk of my PhD project on questions similar to this, and haven't really found a good answer for how to detect unmodelled ancestry.

The general idea of how I tried to approach ancestry determination was to create a model set containing well-defined groups of individuals with a known, specific ancestral background. A query individual (or individuals) were then added to this group, and an ancestry-estimating program was run (specifically Structure) to work out what proportion of the unknown individual was attributable to each of the known groups. I would expect that most genetic ancestry tests follow a similar approach to this, although perhaps with a little bit less care about the definition of the model population groups.

Problems arise when individuals are incorrectly assigned to a particular group, when groups are present in the model set that represent a larger proportion of individuals than most other groups, when groups are closely-related to other groups, and when an ancestral history is present in a test individual that doesn't match any of the model groups. And all that assumes that the marker set used for ancestry determination is perfect: no bias towards any particular group, and no systematic genotyping error.

This doesn't mean that ancestry testing won't work, but it's a good idea to take the results with a large grain of salt. There's a good example of a media personality in New Zealand who was told that she had a very high probability of being 100% Māori, despite having a good knowledge of her own family history that indicated one European ancestor a few generations back on both sides of her family.

# Answer
> 2 votes
Adding another answer that, while much more involved, offers flexibility, walks through minimal steps, and provides visualizations.
I created a Dash app and repository that explores this type of analysis using published ancestry informative SNPs and 1000 Genomes Project Data.

1. Identify published ancestry informative SNPs (AISNPs).
2. Download 1000 Genomes Project genotype bcf data.
3. Limit genotype data from (2) to AISNPs loci from (1).
4. One-hot encode the genotypes.
5. Perform dimensionality reduction (PCA, t-SNE, or UMAP)
6. Plot components.

The original post seemed to want to project their genome against a reference sample. The plot\_walkthrough.ipynb in the `tgviz` repo shows how this can be achieved.
Output from the notebook on my own data:

# Answer
> 1 votes
Related: ethnicity check either from bam or vcf files

The missing piece of your puzzle is an ancestry-labeled reference dataset. 1000 Genomes phase 3 usually works well enough here, and it now has variant callsets for both GRCh37 and GRCh38.

---
Tags: vcf
---


# Question
Title: How do I re-name the headers of my Fasta file?
I appologize for asking this, but I really am really bad with regex...

Can someone help me transform the headers of my fasta files from this:

```
>Agaricus_chiangmaiensis|JF514531|SH174817.07FU|reps|k__Fungi;p__Basidiomycota;c__Agaricomycetes;o__Agaricales;f__Agaricaceae;g__Agaricus;s__Agaricus_chiangmaiensis

```

to this:

```
>JF514531|SH174817.07FU;tax=d:Fungi,p:Basidiomycota,c:Agaricomycetes,o:Agaricales,f:Agaricaceae,g:Agaricus,s:Agaricus_chiangmaiensis_SH174817.07FU;

```

seems like I could use the command

```
sed 's/__/:/g;s/;/,/g;s/.*|/|/g'

```

to substitute `:` for `__`, and `,` for `;`. I'm not entirely sure how to proceed from there though...

Thank you

# Answer
> 5 votes
If `awk` is good as well:

```
$ awk -v FS="|" -v OFS="|" '/>/{
    gsub(";",",", $5);
    gsub("__",":",$5);
    print ">"$2,$3";tax="$5"_"$3";";
    next;
}1' file.fa
>JF514531|SH174817.07FU;tax=k:Fungi,p:Basidiomycota,c:Agaricomycetes,o:Agaricales,f:Agaricaceae,g:Agaricus,s:Agaricus_chiangmaiensis_SH174817.07FU;

```

This relies on the order of the parts separated by `|`.

I don't understand why `k__Fungi` becomes `d:Fungi`, so I ignored it for now.

Some explanation for what `awk`is doing.

* `|` is defined as a field (column) separator for the input and output
* If the line contains a `>` the code between `{` and `}` is executed.
* the first `gsub` replaces `;` with `,` in the 5th field
* the second `gsub` replaces `__` with `:` in the 5th field
* then we print the new header line. The `$2`, `$3` and `$5` get replaced by the values of the second, third and 5th field of the input. `,` get replaced by the defined output field separator
* with `next` we tell `awk` to go to the next input line and repeat the code
* The `1` is in equivalent to `print $0` for printing the whole input line. In our code this is only executed if the line is no header line.

# Answer
> 3 votes
It seems like you want to do four things:

1. Remove all text up to and including the first `|`.
2. Now that that has been removed, take the string between what is now (after step 1) the first `|` and what is now the second `|`, and append it to the end of the header with a `_`.
3. Replace all `__` with `:` and all `;` with `,`.
4. Change the (fixed?) string `reps|k` with `tax=d`.

If so, and assuming that you only have one occurrence of `reps|` in each header, you can do:

```
$ sed 's/>[^|]*|/>/;
       s/reps|k/tax=d/;
       s/__/:/g;
       s/;/,/g;
       s/|\([^|]*\)\(|.*\)/|\1\2_\1;/; ' file.fa
>JF514531|SH174817.07FU;tax=d:Fungi,p:Basidiomycota,c:Agaricomycetes,o:Agaricales,f:Agaricaceae,g:Agaricus,s:Agaricus_chiangmaiensis_SH174817.07FU;

```

This will:

* `s/>[^|]*|/>/` : replace the longest stretch of non`-|` characters after a `>` and until the first `|` with `>` alone.
* `s/reps|k/tax=d/;` : replace the first occurrence of the string `reps|k` with `tax=`.
* `s/__/:/g; s/;/,/g` : replace *all* (that's what the `/g` means) occurrences of `__` with `:` and of `;` with `,`.
* `s/|\([^|]*\)\(|.*\)/|\1\2_\1;/` : the `\(foo\)` construct is called a capturing group and will *capture* `foo` and save it as `\1`, which can be used on the right hand side of the replacement operator. The next captured group will be `\2` then `\3` etc. Here, with `|\([^|]*\)` we are capturing the longest string of non`|` found after the 1st `|` (this happens after the first operation, so the 1st `|` is the one that was 3rd in the original) and then everything after the next `|` (`|.*`) until the end of the line. We then replace everything matched with a `|`, then whatever was matched by the first capture group (`\1`), then the second (`\2`), then we add a `_` and repeat the `\1`. Finally, if you want to make the changes to the original input file, use `-i`:

```
sed -i.bak 's/>[^|]*|/>/; s/reps|k/tax=d/; s/__/:/g; s/;/,/g; s/|\([^|]*\)\(|.*\)/|\1\2_\1;/; ' file.fa

```

That will save the original file as `file.fa.bak`. use `-i` alone without any argument if you don't need the backup.

---
Tags: fasta, text-processing
---


# Question
Title: Parallel in PLINK for linear association for SNP effects
I am doing linear association calculations for SNP effects, I see it take a lot of time to do, does anyone have any experience in parallelizing this job. I have used `--threads <num of cores>` but it says it only uses 1 thread

here is my code for the calculation in plink:

```
plink --bfile ukbb_dis --pheno phenotype_dis.pheno --mpheno 4 --allow-no-sex --covar phenotype_dis.pheno --covar-number 5-21 --missing-phenotype -10000000 --linear --adjust --out assoc_SNP_height_linear --threads 4

```

Here is the current output:

```
257682 MB RAM detected; reserving 128841 MB for main workspace.
1133273 variants loaded from .bim file.
284516 people (0 males, 0 females, 284516 ambiguous) loaded from .fam.
Ambiguous sex IDs written to assoc_SNP_height_linear.nosex .
284516 phenotype values present after --pheno.
Using 1 thread.
Warning: This run includes BLAS/LAPACK linear algebra operations which
currently disregard the --threads limit.  If this is problematic, you may want
to recompile against single-threaded BLAS/LAPACK.
--covar: 17 out of 21 covariates loaded.
Before main variant filters, 284516 founders and 0 nonfounders present.
Calculating allele frequencies... done.
Total genotyping rate is 0.992737.
1133273 variants and 284516 people pass filters and QC.
Phenotype data is quantitative.
Writing linear model association results to
assoc_SNP_height_linear.assoc.linear ... 0%

```

Your help is really appreciated!

# Answer
Try the latest plink 2.0 alpha build for this. It’s an incomplete program, but its implementation of —linear is far better than v1.9’s.
> 2 votes

---
Tags: gwas, plink
---


title: Extracting some part of a list

text:
# Question
Title: Extracting some part of a list
I have two normalized gene expression values (log2 of cpm)

```
> head(biomarker[,1:5])
               A2        A3        A4        A6        A7
A2M    12.4071618 12.631601 12.889748 11.842511 11.574134
ABCB11  0.5151957  6.176788  5.414905  7.134915  1.247590
ABCC2   6.7244303  6.487794 10.178615  6.384132  6.089462
ABCC6   6.3977122  5.823627  7.370761  7.397091  6.071587
ABCF1   9.1609847  9.893258 10.116638 10.225520 10.838486
ABCG2   7.6841874  7.738293  6.833152  6.888041  6.030032
>

```

And

```
> head(immune[,1:5])
              A2        A3        A4        A6        A7
A2M    10.904748 11.388404  9.910614 11.513439 12.963609
ABCB11  5.011380  6.359443  7.145992  8.451947  7.722605
ABCC2   5.040461  6.477014  3.873996  6.409777  9.133971
ABCC6   7.798441  7.601848  9.948072 11.628533 12.701460
ABCF1   8.553597  8.615425 11.145903 10.289098 11.444140
ABCG2   6.224294  5.629375  8.293416  7.979859  8.603793
>

> dim(biomarker)
[1] 719  56
> dim(immune)
[1] 719  56
>

```

I have matched samples in these experiments and same gene too. I could roughly say in both experiments everything has been constant. However, I want to do Student's t-Test to find genes that they had inconsistent gene expression patterns between the two data sets.

By this function I am getting a list but I don't know how to extract the names of genes with p-value \> 0.05 or whatever cut-off

```
> f <- function(x,y){
+     test <- t.test(x,y, paired=TRUE)
+     out <- data.frame(stat = test$statistic,
+                       df   = test$parameter,
+                       pval = test$p.value,
+                       conl = test$conf.int[1],
+                       conh = test$conf.int[2]
+     )
+     return(out)
+ }
> t=sapply(seq(ncol(t(biomarker))), function(x) f(t(biomarker)[,x], t(immune)[,x]))

```

Any help please?

```
  str(t)

    [list output truncated]
     - attr(*, "dim")= int [1:2] 5 719
     - attr(*, "dimnames")=List of 2
      ..$ : chr [1:5] "stat" "df" "pval" "conl" ...
  ..$ : NULL
    >

> sapply(t, "[[", "pval")
Error in FUN(X[[i]], ...) : subscript out of bounds
>

> pvalueGenes <- sapply(colnames(biomarker),function(i) t.test(biomarker[i, ],immune[i, ], paired = TRUE)$p.value)
Error in if (stderr < 10 * .Machine$double.eps * abs(mx)) stop("data are essentially constant") :
  missing value where TRUE/FALSE needed
In addition: Warning messages:
1: In mean.default(x) : argument is not numeric or logical: returning NA
2: In if (stderr < 10 * .Machine$double.eps * abs(mx)) stop("data are essentially constant") :

 Show Traceback

 Rerun with Debug
 Error in if (stderr < 10 * .Machine$double.eps * abs(mx)) stop("data are essentially constant") :
  missing value where TRUE/FALSE needed

```

# Answer
> 2 votes
Starting from the definition of `t`:

```
...
> t <- sapply(seq(ncol(t(biomarker))), function(x) f(t(biomarker)[,x], t(immune)[,x]))

```

Set the column names to the gene names from one or the other of your conditions:

```
> colnames(t) <- rownames(immune)

```

Set the cutoff:

```
> cutoff <- 0.05

```

Filter for rows which meet the p-value cutoff criterium, and get their column names (gene names):

```
> colnames(t)[which(unlist(lapply(t["pval",], function(x) x < cutoff)))]
[1] "ABCC6"
...

```

Only one gene name is reported here, using your demo datasets. Adjust your cutoff and criterium logic, as needed. For example, I filter here for gene names where the p-value is *less* than the cutoff; however, from your question, you may be asking for the inverse, for whatever reason, so you could change `<` to `>=`.

# Answer
> 2 votes
From your example we are comparing samples not genes, I think we need to run the test row wise, something like this:

```
pvalueGenes <- sapply(colnames(biomarker),
                      function(i) t.test(biomarker[i, ],
                                         immune[i, ], paired = TRUE)$p.value)

names(pvalueGenes)[ pvalueGenes < 0.05 ]

```

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


# Question
Title: Albacore basecalling running but outputs 0 reads
I am trying to basecall data produced by the MinION using the SQK-LSK109 kit and FLO-MIN106 flowcell via the command-line. My version of albacore is the latest (v2.3.4). I tried running using the following command-line (the name of the sequencing run and fast5s are different and are default from machine output):

```
python3.6 /home/dkennetz/.local/bin/read_fast5_basecaller.py -i ./seqrun/fast5/ -t 1 -s ./good_reads/ -f FLO-MIN106 -k SQK-LSK109 -r -o fastq

```

-t 1 is just for this test. It is running and finishing almost instantly with 0 errors but outputting no reads. I have previously run this same command-line using the SQK-LSK108 kit and had no problems.

I then went on to validate the fast5s using the h5\_validate tool:

```
python3.6 /home/dkennetz/.local/bin/h5_validate /home/dkennetz/.local/lib/python3.6/site-packages/h5_validator/schemas/single_read_fast5.yaml seqrun/fast5/0/seqrun.fast5 -v

```

The output from this is:

```
    Matcher <'UniqueGlobalKey', type=group match_type=exact count=
{'minimum_count': 1, 'maximum_count': 1}> was not satisfied after matching <HDF5 file "seqrun.fast5" (mode r)>

Error at /:

Matcher <'PreviousReadInfo', type=group match_type=exact count=
{'minimum_count': 1, 'maximum_count': 1}> was not satisfied after matching <HDF5 file "seqrun.fast5" (mode r)>

```

I read on the community page that this may be because the header is missing some relevant meta-data information? Any help on how to fix this would be greatly appreciated!

This is cross-posted on the community page for ONT here:

https://community.nanoporetech.com/posts/albacore-basecalling-runni

# Answer
> 5 votes
I believe you have multiread fast5 files, as generated by recent versions of MinKNOW. A conversion script is available to convert these to the older, one-read one-file format. AFAIK Albacore is being deprecated and no version will be made supporting the multiread fast5 files, as Albacore is being replaced by Guppy.

The command line to convert from multi read fast5s to single read fast5s is as follows:

```
$multi_to_single_fast5 -i /path/to/fast5/ -s /path/to/save/output/

```

---
Tags: nanopore, albacore, base-calling
---


# Question
Title: How to check if indels in VCF files are left or right aligned?
I downloaded a VCF file from dbSNP, and I'm curious if the indels in the file are left-aligned for GRCH37 genome. The documentation doesn't say anything.

How can we tell if a VCF file has left or right aligned indels?

# Answer
One possibility is to run bcftools norm on the file. At the end it will print out a statistic about how many variants were realigned. I did this for dbSNPv151 (GRCh38) only for chromosome 22 like this:

```
$ bcftools view -i 'TYPE="indel"' dbSNP_GRCh38.p7_b151.vcf.gz 22|bcftools norm -Ou -f GRCh38.fa > /dev/null
Lines   total/split/realigned/skipped:  901503/0/490/0

```
> 6 votes

# Answer
By definition in the VCF format specification the POS field is left aligned

> POS int32 t 0-based leftmost coordinate

To confirm that you VCF is following the standard you would have to look at the reference and the match up the entries in the VCF. However if it is from a reliable source like dbSNP you can probably assume they follow the spec
> 5 votes

---
Tags: sequence-alignment, vcf, indel
---


# Question
Title: data storage requirements for RNA-seq and WGS
Supposing upcoming RNA-seq for 3 conditions and 3 replications for each, how much is data storage requirements?

# Answer
> 3 votes
The breakdown would be:

* 1.5-2 GB per sample for fastq files.
* ~3GB per BAM file, depending on the aligner

So ~5GB per BAM file, which is about 50GB total. That can vary by maybe 20%, depending on how compressible the results are (I rearrange our fastq files to make them more compressible, so I'm probably underestimating what others would see).

Obviously using Salmon/Kallisto would halve the storage requirements. Using CRAM would also decrease things, though obviously fewer tools support it.

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


# Question
Title: What does liability mean in GWAS heritability?
I am reading about GWAS in heritability. They usually say to calculate heritability on a liability scale but after searching for this word "liability", I still don't understand clearly what does liability mean? Does anyone have any easy-to-understand explanation regarding this term? Thank you very much!

# Answer
> 5 votes
When you use linear mixed models to estimate heritability you assume that the underlying trait is normally distributed which is called a disease liability scale.

For continuous traits this is not a problem but for binary traits, this becomes an issue because you have a 0/1 value for a phenotype and usually there is a higher proportion of cases in the study sample than the general population prevalence of disease, which leads to an ascertainment bias.

The population prevalence also varies between populations, for example the prevalence of malaria in one continent is vastly different from another.

So to make the heritability estimates comparable and also to correct for the ascertainment bias, the observed scale heritabilities for dichotomous traits are transformed to the liability scale taking population prevalence and the case proportion into account.

Some good papers for reading:

Falconer D.S. The inheritance of liability to certain diseases, estimated from the incidence among relatives. Ann. Hum. Genet. 1965;29:51–71 https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3059431/

---
Tags: gwas
---


# Question
Title: How to export peptide-shaker data from commandline?
I would like to export the data from Protein-Shaker proteomics analysis software and use it in downstream analysis. For example the overview data, numbers of validated peptides per protein, confidence, etc. Ideally I would like to generate this output from commandline without the gui and, if possible to plain txt or csv.

# Answer
You have to start PeptideShaker calling `java` with the `-cp` argument an input project file and provide a comma/space-separated list of integers for the reports you want to export:

```
java -cp PeptideShaker-1.16.37.jar eu.isas.peptideshaker.cmd.ReportCLI -in input.cpsx -reports "0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10"`

```

Here is the list of reports you can export:

```
-reports                   Comma separated list of types of report to export.
                           0: Certificate of Analysis,
                           1: Default Hierarchical Report,
                           2: Default PSM Phosphorylation Report,
                           3: Default PSM Report,
                           4: Default PSM Report with non-validated matches,
                           5: Default Peptide Phosphorylation Report,
                           6: Default Peptide Report,
                           7: Default Peptide Report with non-validated matches,
                           8: Default Protein Phosphorylation Report,
                           9: Default Protein Report,
                           10: Default Protein Report with non-validated matches,
                           11: Extended PSM Report,
                           12-n: Your own custom reports.

```

I don't know whether you can also use the SearchGUI output zip to generate this files. The program accepts zip files as input, but I think that means PeptideShaker project that are exported to zip files, but not SearchGUI output zips. So, in the end you might need to run the GUI at least one time, unfortunately.
> 1 votes

---
Tags: peptide-shaker
---


# Question
Title: What is the actual cause of excessive zeroes in single cell RNA-seq data? Is it PCR?
First, sorry if I am missing something basic - I am a programmer recently turned bioinformatician so I still don't know a lot of stuff. This is a cross post with a Biostars question hope that's not bad form.

---

While it is obvious that scRNA-seq data contain lots of zeroes, I couldn't find any detailed explanation of why they occur except for short notices along the lines of "substantial technical and biological noise". For the following text, let's assume we are looking at a single gene that is not differentially expressed across cells.

If zeroes were caused solely by low capture efficiency and sequencing depth, all observed zeroes should be explained by low mean expression across cells. This however does not seem to be the case as the distribution of gene counts across cells often has more zeroes than would be expected from a negative binomial model. For Example the ZIFA paper explicitly uses a zero-inflated negative binomial distribution to model scRNA-seq data. Modelling scRNA-seq as zero-inflated negative binomial seems widespread throughout the literature.

However assuming negative binomial distribution for the original counts (as measured in bulk RNA-seq) and assuming that every RNA fragment of the same gene from every cell has approximately the same (low) chance of being captured and sequenced, the distribution across single cells should still be negative binomial (see this question for related math).

So the only remaining possible cause is that inflated zero counts are caused by PCR. Only non-zero counts (after capture) are amplified and then sequenced, shifting the mean of the observed gene counts away from zero while the pre-PCR zero counts stay zero. Indeed some quick simulations show that such a procedure could occasionally generate zero-inflated negative binomial distributions. This would suggest that excessive zeroes should not be present when UMIs are used - I checked one scRNA-seq dataset with UMIs and it seems to be fit well by plain negative binomial.

Is my reasoning correct? Thanks for any pointers.

The question How can we distinguish between true zero and dropout-zero counts in single-cell RNA-seq? is related, but provides no clues to my present inquiry.

# Answer
It may be necessary to distinguish between methods that use unique molecular identifiers (UMIs), such as 10X's Chromium, Drop-seq, etc, and non-UMI methods, such as SMRT-seq. At least for UMI-based methods, the alternative perspective, that **there is no significant zero-inflation in scRNA-seq**, is also advocated in the single-cell research community. The argument is straight-forward: the empirical mean expression vs. dropout rate curve matches the theoretically predicted one, *given the current levels of capture efficiency*.

# Examples

### Svensson Blog

A couple of blog posts from Valentine Svensson argue this point rather pedagogically, and include citations from across the literature:

Droplet scRNA-seq is not zero-inflated

Count-depth variation makes Poisson scRNA-seq data negative binomial

### bayNorm

There is a more extensive preprint by Tang, Shahrezaei, et al. (BioRxiv, 2018) that claims to show a binomial model is sufficient to account for the observed dropout noise. Here is a snippet of a relevant conclusion:

> Importantly, as bayNorm recovered dropout rates successfully in both UMI-based and non-UMI protocols without the need for specific assumptions, we conclude that invoking zero-inflation models is not required to describe scRNA-seq data. Consistent with this, the differences in mean expression levels of lowly expressed genes observed between bulk and scRNA-seq data, which were suggested to be indicative of zero-inflation, were recovered by our simulated data using the binomial model only.

### Multinomial Modeling

There is also a very clearly written preprint by Townes, Irizarry, et al. (BioRxiv, 2019) where the authors consider scRNA-seq as a proper compositional sampling (i.e., multinomial process) and they come to a similar conclusion, though *specifically for UMI-based methods*. From the paper:

> The multinomial model makes two predictions which we verified using negative control data. First, the fraction of zeros in a sample (cell or droplet) is inversely related to the total number of UMIs in that sample. Second, the probability of an endogenous gene or ERCC spike-in having zero counts is a decreasing function of its mean expression (equations provided in Methods). Both of these predictions were validated by the negative control data (Figure 1). In particular, the empirical probability of a gene being zero across droplets was well calibrated to the theoretical prediction based on the multinomial model. This also demonstrates that UMI counts are not zero inflated.

Furthermore, by comparing raw read counts (prior to UMI-based deduplication) and UMI counts, they conclude that PCR is indeed the cause of zero-inflation:

> The results suggest that while read counts appear zero-inflated and multimodal, UMI counts follow a discrete distribution with no zero inflation (Figure S1). The apparent zero inflation in read counts is a result of PCR duplicates.

I highly recommend giving this a read, especially because it nicely situates other common generative models (e.g., binomial, Poisson) as valid simplifying assumptions of the multinomial model.

It should be noted that this same group previously published a work (Hicks, Irizarry, et al. 2018), mostly focused on non-UMI-based datasets (SMRT-seq), where they showed evidence that, relative to bulk RNA-seq, there was significant zero-inflation.
> 9 votes

# Answer
The Biostars thread turned out helpful. The most interesting possible cause, not mentioned in the Ian Subery's answer, is that due to bursty nature of transcription, the true distribution of transcript counts across cells can be bimodal with a peak at zero even assuming a simple model of transcription such as the random telegraph model. See for example Dattani & Barahona 2017 for a more detailed discussion.

However, the random telegraph model predicts that the even with peak at zero, there should usually be non-negligible probability of having small, but nonzero counts, which is not always the case in zero-inflated scRNA-seq data: see for example the Isl1 gene in Retina dataset (click Explore -\> search for "Isl1")
> 8 votes

# Answer
I know of no references for this, but in general, I would say that your reasoning is sound. I would just add that in contrast to what I suspect you have simulated, not all transcripts are equally likely to be captured and amplified. We don't really understand what the determinants of this are, but for example, GC content is definitely related.
> 5 votes

---
Tags: rna-seq, scrnaseq, single-cell, modelling
---


# Question
Title: How to quantile normalization on RNA seq counts
I have a read count data (RNAseq) and want to perform quantile normalization. Could you please help me how to do it. I tried some scripts in R but it didn't work. I want the result output in matrix form.

```
  gene_id   SRR896664   SRR896663   SRR896665
  ENSG00000000003   46106   36353   40614
  ENSG00000000005   198 399 1200
  ENSG00000000419   40364   37769   40849
  ENSG00000000457   18924   16211   16057
  ENSG00000000460   31040   28888   29901
  ENSG00000000938   200 0   394
  ENSG00000000971   14935   14353   12522

```

The script which i tried is

```
  data <- read.csv("testquantile.csv",header=T)
  head(data)
  rownames(data) <- data[,1]
  data_mat <- data.matrix(data[,-1])
  head(data_mat)
  data_norm <- normalize.quantiles(m, copy = TRUE)

```

# Answer
On google there are many tutorials about quantile normalzation, for example here. In that tutorial they made a function to calculate quantile normalization. Here an example with that function on your small data set.

```
data
          gene_id SRR896664 SRR896663 SRR896665
1 ENSG00000000003     46106     36353     40614
2 ENSG00000000005       198       399      1200
3 ENSG00000000419     40364     37769     40849
4 ENSG00000000457     18924     16211     16057
5 ENSG00000000460     31040     28888     29901
6 ENSG00000000938       200         0       394
7 ENSG00000000971     14935     14353     12522

rownames(data) <- data$gene_id

quantile_normalisation <- function(df){
  df_rank <- apply(df,2,rank,ties.method="min")
  df_sorted <- data.frame(apply(df, 2, sort))
  df_mean <- apply(df_sorted, 1, mean)

  index_to_mean <- function(my_index, my_mean){
    return(my_mean[my_index])
  }

  df_final <- apply(df_rank, 2, index_to_mean, my_mean=df_mean)
  rownames(df_final) <- rownames(df)
  return(df_final)
}

new_data <- quantile_normalisation(data[,2:4])
new_data
                 SRR896664  SRR896663  SRR896665
ENSG00000000003 41574.6667 39110.3333 39110.3333
ENSG00000000005   197.3333   599.6667   599.6667
ENSG00000000419 39110.3333 41574.6667 41574.6667
ENSG00000000457 17064.0000 17064.0000 17064.0000
ENSG00000000460 29943.0000 29943.0000 29943.0000
ENSG00000000938   599.6667   197.3333   197.3333
ENSG00000000971 13936.6667 13936.6667 13936.6667

boxplot(data[,2:4])

```

```
boxplot(new_data)

```
> 4 votes

# Answer
Ma be CQN from Bioconductor will be useful, though it doesn't perform just quantile normalisation.
> 3 votes

# Answer
> https://www.rdocumentation.org/packages/limma/versions/3.28.14/topics/normalizeQuantiles

is a possibility.

Another popular method is https://www.rdocumentation.org/packages/preprocessCore/versions/1.34.0/topics/normalize.quantiles.

```
library(preprocessCore)
normalize.quantiles

```
> 1 votes

---
Tags: rna-seq, gene, normalization, gene-expression
---


# Question
Title: Can a gff file be converted to a fasta file?
I downloaded an annotated genome file in gff format here. I would like to use it for proteomics. Though I need it in fasta format. Is there any tool that converts gff to annotated fasta? I see the file contains the genome sequence in fasta format, but not annotated. Are those file formats convertable?

# Answer
A gff is a file of annotation. It generally doesn't include sequence information, so you can not inter-convert. But the file on your page has a fasta entry at the bottom.. I'd just grab the bottom half of it, and if the sequence length looks right, that's likely what you want.
> 5 votes

# Answer
You can use GFF utilities to get a fasta sequence from a GFF file using a reference genome file.
> 5 votes

# Answer
> Is there any tool that converts gff to annotated fasta?

This is not an uncommon question / request, but perhaps a crash course refresher on file formats would provide some valuable context.

* **FASTA files** are designed to hold nucleotide or amino acid sequences. The "defline" (header) of each record may include some unique identifiers or descriptive metadata. There are no strict requirements, although there are many common conventions. The information that *should* be the defline of a FASTA record depends primarily on the tool(s) used to produce and analyze the data.
* **GFF files** are designed to hold annotations corresponding to sequences, most often genomic sequences such as chromosomes or (if it's a draft genome assembly) scaffolds or contigs. Each record in the GFF file specifies a sequence (such as `chr6` or `scaffold1985`) and an interval (such as `19920557-19997344` or `2488-2489`). The GFF file labels that interval as a gene or transposable element or transcription factor binding site or whatever genomic feature you're interested in.

So in this context, a "*tool that converts gff to annotated fasta*" doesn't make sense since it's an underdefined task. There are many ways to extract sequences from a Fasta file that has been annotated with a GFF file. Many, many scientists have written throwaway programs to do this type of task. It can be done in a hundred subtly different ways, and most of these programs are designed for a single narrow use case. Do you want to:

* pull out the sequence corresponding to every GFF record? Usually not helpful, since many records are children of other records.
* pull out the sequence corresponding to every gene? More common, but since you want to do proteomics you're probably not interested in promoters or introns (if either of these are annotated in your GFF file).
* pull out the coding sequence of each gene and translate it into its corresponding peptide sequence?

It looks like that last point is most likely, but unless I missed something it looks like Prokka is already doing that for you. So the question remains: **precisely what information exists in your GFF file that you would like in your FASTA file**? You've provided an example of what you have (the FASTA files at least), but can you provide an example of what you *want*? Only then can we help you to identify the right tool for the task. There's a good chance it will involve writing some code in a programming language like Python or Perl.
> 3 votes

# Answer
I extracted the sequence contained in the gff file into a separate file and added a header according to the fasta file format specifications. I named that file `ERS654933.fasta` and then used `prokka` to generate my own annotation.

File header ERS654933.fasta:

```
>ATCC_NCTN13373
GTGTATAACTTAAAAATTTAAGAAAGATGGAGTAAATTTATGTCGGAAAAAGAAATTTGG
GAAAAAGTGCTTGAAATTGCTCAAGAAAAATTATCAGCTGTAAGTTACTCAACTTTCCTA
AAAGATACTGAGCTTTACACGATCAAAGATGGTGAAGCTATCGTATTATCGAGTATTCCT
TTTAATGCAAATTGGTTAAATCAACAATATGCTGAAATTATCCAAGCAATCTTATTTGAT
GTTGTAGGCTATGAAGTAAAACCTCACTTTATTACTACTGAAGAATTAGCAAATTATAGT

```

Then:

```
prokka ERS654933.fasta

```

Which generated a folder with several files. Here is the header of the generated file `PROKKA_01232019.faa`.

```
NEOBMGJD_00001 Chromosomal replication initiator protein DnaA
MSEKEIWEKVLEIAQEKLSAVSYSTFLKDTELYTIKDGEAIVLSSIPFNANWLNQQYAEI
IQAILFDVVGYEVKPHFITTEELANYSNNETATPKEATKPSTETTEDNHVLGREQFNAHN
TFDTFVIGPGNRFPHAASLAVAEAPAKAYNPLFIYGGVGLGKTHLMHAIGHHVLDNNPDA
KVIYTSSEKFTNEFIKSIRDNEGEAFRERYRNIDVLLIDDIQFIQNKVQTQEEFFYTFNE
LHQNNKQIVISSDRPPKEIAQLEDRLRSRFEWGLIVDITPPDYETRMAILQKKIEEEKLD
IPPEALNYIANQIQSNIRELEGALTRLLAYSQLLGKPITTELTAEALKDIIQAPKSKKIT
IQDIQKIVGQYYNVRIEDFSAKKRTKSIAYPRQIAMYLSRELTDFSLPKIGEEFGGRDHT
TVIHAHEKISKDLKEDPIFKQEVENLEKEIRNV
>NEOBMGJD_00002 Beta sliding clamp
MMEFTIKRDYFITQLNDTLKAISPRTTLPILTGIKIDAKEHEVILTGSDSEISIEITIPK
TVDGEDIVNISETGSVVLPGRFFVDIIKKLPGKDVKLSTNEQFQTLITSGHSEFNLSGLD
PDQYPLLPQVSRDDAIQLSVKVLKNVIAQTNFAVSTSETRPVLTGVNWLIQENELICTAT
DSHRLAVRKLQLEDVSENKNVIIPGKALAELNKIMSDNEEDIDIFFASNQVLFKVGNVNF
ISRLLEGHYPDTTRLFPENYEIKLSIDNGEFYHAIDRASLLAREGGNNVIKLSTGDDVVE
LSSTSPEIGTVKEEVDANDVEGGSLKISFNSKYMMDALKAIDNDEVEVEFFGTMKPFILK
PKGDDSVTQLILPIRTY
>NEOBMGJD_00003 hypothetical protein
MVQEVVVEGDINLGQFLKTEGIIESGGQAKWFLQDVEVLINGVRETRRGKKLEHQDRIDI
PELPEDAGSFLIIHQGEQ

```

So, that way I could create my own annotation, but I am loosing the annotation information in the `gff` file.
> 2 votes

---
Tags: fasta, file-formats, format-conversion, gff
---


title: scRNA-seq differential transcript usage

text:
# Question
Title: scRNA-seq differential transcript usage
Many of the modern gene-quantification tools (Salmon/Kallisto) output transcript-level (as opposed to gene-level) data. All of the scRNA-seq analysis I have seen just uses the gene-level values. I suppose I can just use the transcript-level matrix, but then I lose any links between the different isoforms of the same gene. Maybe it does not really matter?

There are a few options for bulk RNA-seq and there was a similar question here, but are there any single-cell options?

# Answer
> 4 votes
> I was hoping to find some tool that already exists.

The short answer is, **none exist yet**. I see a number of hurdles to overcome first.1

For example, **no one has yet characterized saturation for transcript discovery**, i.e., the reads per cell needed to maximize transcript information recovery. This is a typical stepping-stone as a sequencing approach matures, and until it's done, it makes it unlikely that most scRNA-seq datasets are going to be appropriate vis-à-vis statistical power, even if there was a tool. That is, if people are only aiming for cell-type discrimination via gene-level markers, datasets will likely be underpowered to measure transcript level differences.

Another issue is the **lack of any transcript-level scRNA-seq simulators**. All the simulators out there that I know of produce gene-level matrices. Having a transcript-level simulator is crucial for benchmarking and characterizing the performance of differential expression calling methodologies. Without anything to test against, that leaves a huge blindspot in terms of estimating the uncertainty (false positive rates, sensitivity) that any tool might have. However, since some scRNA-seq models (e.g., Velten, et al. 2015, or Lopez, et al., 2018) are generative, perhaps this is shouldn't be too difficult to adapt to generating FASTQs.

Simulators are fine and all, but ground truth is probably the key hurdle of the enterprise right now. While there are definitely groups using smFISH to distinguish isoforms quantitatively on a single-cell level, **there doesn't appear to be enough smFISH data out there to measure a transcriptome-wide scRNA-seq method against**. Hopefully this will change in the coming year or two, as results from more highly-multiplexed single-molecule technologies become available.

**1** This is just my opinion formed as I try to keep up with the literature, so I could be missing things. Comments/corrections certainly welcome!

# Answer
> 2 votes
There is a paper that takes advantage of transcript compatibility counts (TCCs). A TCC describes all of the potential transcripts that a cDNA alignment can be derived from. The senior author, Lior Pachter also has a nice blog post about differential transcript usage.

A discriminative learning approach to differential expression analysis for single-cell RNA-seq

> Single-cell RNA-seq makes it possible to characterize the transcriptomes of cell types across different conditions and to identify their transcriptional signatures via differential analysis. Our method detects changes in transcript dynamics and in overall gene abundance in large numbers of cells to determine differential expression. When applied to transcript compatibility counts obtained via pseudoalignment, our approach provides a quantification-free analysis of 3′ single-cell RNA-seq that can identify previously undetectable marker genes.

---
Tags: scrnaseq, single-cell, transcriptome
---