# Question
Title: Merge 2 VCFs from different variant callers
I am working with WES data for detection of somatic variants and I have used two variant callers because no variant caller is complete in itself. I have used GATK Haplotypecaller for small variants like snp,ins,del etc. And for large variants I have used pindel (specially for SVs detection). Now I have two vcf file for each patients, one from gatk and another from pindel. Although I know various vcf-merging tools but is it right to merge two different vcf file from different callers. Do they give right and proper merged vcf file. Are those vcf-merging tools SV-aware and consider only trivial call for the same variant represented differently in different callers. I know these are too many questions for a single thread but all these are related to each other. I am in dilemma whether I should merge them or use them separately.

Thanks.

# Answer
The merge utility from SURVIVOR seems to be what you are looking for. I did not try it yet, but it seems to be especially designed to handle SV.

The simplest thing to do would probably to first merge your patients files together for each tool and have one multi-patient file per variant caller. Then I would merge the multi-sample files of the different tools.

Here is a detailed blog post describing the different challenges associated with merging structural variants from different callers and explaining the steps they took to address it. All the code to reproduce the pipeline is available on a github repo.
> 4 votes

---
Tags: vcf, variant-calling, gatk, variation
---


# Question
Title: Get gene names corresponding to CPG probes
I have downloaded KIRC Methylation data (450K) using TCGA toolbox. The beta values given are with reference to CPG probes. Is there any package or function in R through which I can get the gene names corresponding to these probes? I was able to get the gene I.Ds using `getMappedEntrezIDs`. Is there any function to map these I.Ds to respective gene names?

# Answer
Here's a Python script (e.g., called `map_entrez_to_ensembl.py` that uses mygene.info to convert a list of Entrez IDs to Ensembl gene names:

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

import sys
import mygene

mg = mygene.MyGeneInfo()

genes = ['3815',
         '3816',
         '2341']

results = mg.querymany(genes, scopes=['entrezgene'], fields=['ensembl'], species='human', verbose=False)

for res in results:
    entrez_id = res['query']
    if 'ensembl' in res and 'gene' in res['ensembl']:
        sys.stdout.write('%s\t%s\n' % (entrez_id, res['ensembl']['gene']))

```

Result:

```
$ ./map_entrez_to_ensembl.py
3815    ENSG00000157404
3816    ENSG00000167748

```

Looks like there's an R-based interface to mygene.info, if you really need to use R.

You could use their API test page to see how data are structured and what keywords are available for different kinds of searches.

For instance, if you wanted aliases, in addition to Ensembl names:

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

import sys
import mygene

mg = mygene.MyGeneInfo()

genes = ['3815',
         '3816',
         '2341']

results = mg.querymany(genes, scopes=['entrezgene'], fields=['ensembl', 'alias'], species='human', verbose=False)

for res in results:
    entrez_id = res['query']
    if 'ensembl' in res and 'gene' in res['ensembl'] and 'alias' in res:
        for alias in res['alias']:
            sys.stdout.write('%s\t%s\t%s\n' % (entrez_id, alias, res['ensembl']['gene']))

```

Then:

```
$ ./map_entrez_to_alias_and_ensembl.py
3815    C-Kit   ENSG00000157404
3815    CD117   ENSG00000157404
3815    PBT     ENSG00000157404
3815    SCFR    ENSG00000157404
3816    KLKR    ENSG00000167748
3816    Klk6    ENSG00000167748
3816    hK1     ENSG00000167748

```
> 3 votes

---
Tags: r, illumina
---


# Question
Title: t-sne code and labelling issue for rna seq data
I'm trying to run t-sne on my coding as well as non coding genes. Even though I'm getting output I'm not confident about my code (The output it's difficult to interpret).

I have a data set of total 16 cell type and about 3000 genes.

This is how my data is first column is gene rest is sample:

```
gene    H1  H2  H3  H4  M2  M3  M4  C1  C2  C3  C4  G2  G3  G4  Gran1   Gran2

```

The code used:

```
library(Rtsne)
file1<- read.csv('UNIQUE_WGCNA_FILE.csv',header = T,row.names = 1)
trn <- data.matrix(file1)
dim(trn)

tsne <- Rtsne(as.matrix(trn[,1:16]), check_duplicates = FALSE, pca = FALSE, perplexity=50, theta=0.5, dims=3)

cols <- rainbow(10)
plot(tsne$Y, t='n')
text(tsne$Y, labels=as.numeric(trn[,1]), col=cols[trn[,1]])

require(rgl)
plot3d(tsne$Y, col=cols[trn[,1]])
legend3d("topright", legend = '0':'5', pch = 16, col = rainbow(10))
rgl.postscript("Merge.pdf", "pdf")

```

I'm getting output but I don't know how to label so its becoming difficult. Any suggestion or help how to label would be appreciated.

# Answer
> 3 votes
Although a bit late, here is the solution for this:

1) Create color/column levels:

```
levels = unlist(lapply(colnames(trn[,1:16]), function(x)
    {

        if(x == "H1") return("H1")
        else if(x == "H2") return("H2")

}))

```

2)

```
plot(tsne$Y, col=factor(levels))

```

This should color columns H1 and H2 differently

Edit:

This is a sample complete code:

```
library(Rtsne)
file1<- read.csv('filename.csv',header = T, row.names=1)
trn <- data.matrix(file1)
dim(trn)

tsne <- Rtsne(as.matrix(trn[,1:16]), check_duplicates = FALSE, pca = FALSE, perplexity=5, theta=0.5, dims=3)

collevels = unlist(lapply(colnames(trn), function(x)
{
  if(grepl("H1",x))
    return("blue")
  else if(grepl("H2",x))
    return("red")

  # etc ... for all the columns you want to color
}
))

colorColumns = factor(collevels)
require(rgl)
plot3d(tsne$Y[,1], tsne$Y[,2], tsne$Y[,3], col=colorColumns)

```

# Answer
> 1 votes
It seems that `tsne$Y` is the matrix and the names are just taken from there. So, you could replace the names of the columns of that matrix and then it would display the names instead of the numbers.

But, you are using the function `plot3d` from the rgl package, which has as arguments `xlab`, `ylab`, and `zlab`, so you just need to give a string to those arguments to name the axes.

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


title: How to add data to AnnotationHub

text:
# Question
Title: How to add data to AnnotationHub
I recently started using the Bioconductor AnnotationHub package for accessing lots of common datasets, but I've noticed that there are a lot of newer databases that have not made it into AnnotationHub's collection.

Has anyone had any luck getting databases added to AnnotationHub or adding them manually yourself (I guess via the package process?)? Is it painful or relatively straightforward?

I love the annotationhub interface, but I'm wondering if it is way too inefficient to add new resources rather than just manually downloading and curating myself. Thanks!

# Answer
As the description of the package says:

> The AnnotationHub web resource provides a central location where genomic files (e.g., VCF, bed, wig) and other resources from standard locations (e.g., UCSC, Ensembl) can be discovered.

If the new resources aren't in the standard locations it won't be listed by AnnotationHub. The easy/common way to have the data accessible through AnnotationHub is make it available through these standard locations (before publications if you want to give early access).

Another method would be to follow this steps of the vignette to add new resources, but it might involve creating a new package.
> 2 votes

---
Tags: bioconductor, sequence-annotation, public-databases
---


# Question
Title: Will blastn remove sequences from a search with low identity?
I'm using this command (excuse the duplicate naming, I know it's bad form):

```
blastn -query cm-seqs/combined_seqs.fna -db combined_seqs.fna -out cm-matched.txt -num_alignments 1 -outfmt 10

```

to take a set of fasta sequences to blast against another data base I built with blast.

My query file has 4,364,417 sequences. The files resulting from the above command spat out 4,362,639 results. The difference is about 1800. Is this due to poor alignments for some sequences that didn't pass any of blasts default parameters? I assumed the result would be 1 "best" match in the database for each query sequence.

# Answer
> 1 votes
Yes, if a query sequence has no good alignment, blast will not return a result for it. Your command will indeed output only the best hit for each query sequence, but there is no reason to assume *all* of your sequences will have a match.

Those that don't have a match won't be included in the output. This sounds perfectly normal unless you are blasting the same input file against a DB created from that file. Otherwise, you don't expect all of your sequences to necessarily have a match.

Also note that since the default blastn settings are quite permissive (usually an e-value of 10, for example), don't make the mistake of assuming that a blast hit means anything. The details will always depend on the specific search you are running (I have discussed this a little in my answer here), but in the vast majority of cases, any hit with an e-value \> 1 can safely be ignored. In most cases, you are only interested in e-values well under 1, actually. So you probably already have all sorts of spurious results in there anyway.

For more details, please edit your question and explain what the query and DB files are.

---
Tags: blast
---


title: Human SNP databases other than dbSNP

text:
# Question
Title: Human SNP databases other than dbSNP
I'm working on a project to identify variants as being either germline or somatic based on nearby SNPs, and I was wondering if there were any databases of SNPs for Homo sapiens other than dbSNP? I'm not looking for databases that contain additional information (pathogenicity etc.) on SNPs that are already in dbSNP, instead I want to know if there are any other databases that contain SNPs that aren't found in dbSNP.

I've been searching online and so far the only databases I can find are either deprecated or integrated into dbSNP, so any suggestions would be great.

# Answer
I'm pretty sure there are no public databases of variants that aren't included in dbSNP. That is *the* variant database so that's the one you want to use.

While there are various other databases around, they tend to use dbSNP rsIDs as well so any variants in those will also be in dbSNP. So while dbSNP is nowhere near complete, of course, it's still the most complete one you're likely to find.
> 3 votes

---
Tags: database, public-databases, snp, variant-calling
---


# Question
Title: Convert CD (cluster of differentiation) into gene symbol
I want to convert names according to the "cluster of differentiation" nomenclature into accepted gene symbols. For example, CD3 should map to "CD3a;CD3b;CD3c;CD247" and "CD45RO" should map to "PTPRC". If the surface antigen isn't coded for by a gene, it should map to "" or something similar.

I've seen Obtain CD (cluster of differentiation) marker expression information, but the Abcam table used there does not give the gene symbol -- it gives a list of alternate names. I've seen table 1 of this paper, but it contains errors -- for example, TNFRSF11 should be TNFRSF11A. I could copy this table by hand, or I could download this one, though that's annoying because it doesn't have one row per CD.

Is there a canonical (or just a convenient and reliable) way to do this?

# Answer
> 2 votes
I ended up using option 4 and parsing the columns for "Synonym" and "Previous Symbol".

https://www.genenames.org/cgi-bin/genefamilies/set/471

---
Tags: public-databases
---


# Question
Title: Declaring an array of strings and integers
I have two csv file - one is a list of strings and other has both strings and floating values. If the strings in both the files match I need to copy the complete row from the second file to a temporary array. I have declared a matrix with n rows and m columns and saved is as dataframe `temp <- as.data.frame(matrix(nrow = n, ncol = m))`. I am using `for` loop to iterate over all rows and `strcmp` to compare the strings. But when the data is written to `temp` the strings are stored as integers. How should I declare this array? Is it possible to write this loop without declaring an array in R?

# Answer
The `merge` function of R might be useful. The function expects two data frames, with \[optionally\] columns to compare, and it will output all fields from both data frames where the comparison columns match.

```
> a <- data.frame(test=1:10);
> b <- data.frame(test=c(sample(1:10,5),sample(c("A","B","C","D","E"))),
                  output=rnorm(10));
> merge(a, b, by="test");
  test     output
1    1 -1.3800443
2    2  0.2934274
3    3  0.0868700
4    9  1.1766914
5   10  0.3793448

```

The `merge` function can be customised quite a bit (e.g. if you want NA set for non-matching rows), see `?merge` for more details.
> 1 votes

# Answer
When you read in your two CSV files, you can explicitly declare the variable types for each column:

```
> a <- read.csv(..., colClasses=c("character"), ...)
> b <- read.csv(..., colClasses=c("character", "numeric", "numeric", ...), ...)

```

You could also investigate `as.is` or `stringsAsFactors`, although this can cause problems where there are duplicate strings. Using `colClasses` is more explicit.

Read the help documentation in `?read.csv` for a more complete explanation of these operators.

You can filter with the `%in%` operator in `match`:

```
> strings_of_interest <- as.vector(a$V1)
> b.filtered <- b[b$V1 %in% strings_of_interest,]

```

If I am reading your question correctly, perhaps you can then take `b.filtered` and turn that into the matrix you want.

Avoid use of `for` loops where you can. They can slow R down considerably. A good read on this subject is in the third chapter of Patrick Burns' The R Inferno.
> 1 votes

---
Tags: r
---


# Question
Title: Keep Format and Individual fields when annotating VCF with VEP
I'm currently updating my Variant Calling Pipeline by switching the VCF annotating software from Annovar to VEP for a variety of regions, not least how easy it is to annotate with HGVS notation and keep datasets up to date in VEP.

For the most part everything is running smoothly, with the exception that some of the data in the VCFs is lost during annotation (and conversion to tsv). The VCFs are created with GATK's UnifiedGenotyper and include a 'Format' column where each value is 'GT:AD:DP:GQ:PL' and a column named after the Individual, which contains semicolon-separated data that corresponds to the Format column (i.e. Genotype;Allele Depth;Depth;Genotype Quality;Phred-likelihood). When I annotate with VEP none of this data is carried over to the output file as it would be in Annovar, leaving me with an annotated file that has no information on read depth, genotype or any of the other data in the two lost columns.

I've included the command I'm currently using for annotation:

```
./vep -i RM0108.vcf --cache --force_overwrite --tab --merged --variant_class --sift b --polyphen b --hgvs --symbol --canonical --check_existing --af_1kg --af_gnomad --humdiv --pick -o RM0108.tsv

```

I can't find information on this in the VEP documentation or elsewhere online. I could write something to take the relevant information from the VCF and add it to the tsv after VEP has finished running, but it seems like there may be an easier solution that I'm missing, so any help would be appreciated.

I've also posted this question on BioStars Link Here

# Answer
The VEP TSV format can only keep in its own specified columns. If you want to maintain the data from your original input, get your output in VCF. It will add the VEP annotation to the INFO column, and keep all the stuff you already have there.
> 4 votes

---
Tags: vcf, sequence-annotation, variant-calling
---


title: What does the setSize ordering mean?

text:
# Question
Title: What does the setSize ordering mean?
I am running `gseGO` function from `clusterProfiler` and getting an output where terms are ordered by the `setSize` parameter. What does this ranking mean? Are terms ordered by the sizes of the gene sets? If yes, then what is the ultimate value of such ranking? My goal is to get the most relevant terms, but I am not completely sure that such an ordering based on the `setSize` is actually the right one.

# Answer
As you can see on this blog post of the author, the result is ordered by `pvalue`, but it can coincidentally sort by `setSize`. The `setSize` means that that term has so many genes on under that category.
> 0 votes

---
Tags: r, clusterprofiler
---


title: Finding proteins in a sequence

text:
# Question
Title: Finding proteins in a sequence
I have a file with protein sequences, all on a single line, like this:

```
>Contig_1 1FFPMPQ*CCLFSSSYLPA*AIMVTGSCQTLVQDKV*CLASCHVMIPPCVFHPIYPLFTWSPLHPIISGCF*AVCMCSWIHALFVEIILSWSLSLPFTFEISSD>Contig_2 1FFPCHNSVVYFLLLICLLKPLWSQGVVRHWYRTRCDAWLAVML*YLHVSFIPSIHFLLGRLCTLLYLDVFELYACAAGFMPCLLR*FSLGHSLFLLPLKFPATVKFFSLLGQIMNNISLMNDLQTGEKLQETH*EVEDKSALLVSIIESRCCVTLQPQDFSN*HLGKL*KSVGQKLK*VNDQITY*

```

I want to print all sequences between an `M` residue and a `*` that are at least 6 residues long, and I would like the output like this:

```
>contig_1.1.1
M.......*
>contig_1.1.2
M........*
>contig_2.1.1
M............*
>contig_2.1.2
M......*

```

Where the first number after the underscore is the number found in as the first character of the sequence (it is the reading frame) and the second number is the sequence number. and the last number is the fragment number. So, since the 1st sequence starts with;

```
>Contig_1 1FFPMPQ

```

The name will be `>Contig_1.1` and each fragment found will be `>Contig_1.1.1`, `>Contig_1.1.2`, ... , `>Contig_1.1.N`.

After implementing the suggestion in the answer below and escaping the `*` in my regular expression, my current code is:

```
import sys
import re

myprotein = ">Contig_1 1FFPMPQ*CCLFSSSYLPA*AIMVTGSCQTLVQDKV*CLASCHVMIPPCVFHPIYPLFTWSPLHPIISGCF*AVCMCSWIHALFVEIILSWSLSLPFTFEISSD>Contig_2 1FFPCHNSVVYFLLLICLLKPLWSQGVVRHWYRTRCDAWLAVML*YLHVSFIPSIHFLLGRLCTLLYLDVFELYACAAGFMPCLLR*FSLGHSLFLLPLKFPATVKFFSLLGQIMNNISLMNDLQTGEKLQETH*EVEDKSALLVSIIESRCCVTLQPQDFSN*HLGKL*KSVGQKLK*VNDQITY*"
# parse record
my_protein = []
for line in myprotein:
    if not line.startswith('>'):
        my_protein.append(line)
    else:
        header = line

# build sequence
my_pro = ''.join(my_protein)

# cleave sequence
peptides = re.findall(r'(M[A-Z]{6,}\*)', my_pro)

# build a formatted list
records = ["%sNo.%d\n%s\n" % (header, idx + 1, fragment) for idx, fragment in     enumerate(peptides)]

# print list
for record in records:
    sys.stdout.write("%s" % (record))

```

And it prints:

```
 >No.1
 MVTGSCQTLVQDKV*
 >No.2
 MIPPCVFHPIYPLFTWSPLHPIISGCF*
 >No.3
 MNNISLMNDLQTGEKLQETH*

```

How can I change this to get my desired output format?

# Answer
> 2 votes
The `*` character has a special meaning in regular expressions, it means "0 or more". So when you write `{6,}*`, the regex engine tries to apply the "0 or more" to the `{6,*}` and gets confused:

```
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python3.6/re.py", line 222, in findall
    return _compile(pattern, flags).findall(string)
  File "/usr/lib/python3.6/re.py", line 301, in _compile
    p = sre_compile.compile(pattern, flags)
  File "/usr/lib/python3.6/sre_compile.py", line 562, in compile
    p = sre_parse.parse(p, flags)
  File "/usr/lib/python3.6/sre_parse.py", line 855, in parse
    p = _parse_sub(source, pattern, flags & SRE_FLAG_VERBOSE, 0)
  File "/usr/lib/python3.6/sre_parse.py", line 416, in _parse_sub
    not nested and not items))
  File "/usr/lib/python3.6/sre_parse.py", line 619, in _parse
    source.tell() - here + len(this))
sre_constants.error: multiple repeat at position 10

```

In order to use special characters literally, as simple strings, you need to escape them by adding a `\` before them. So try:

```
import re
myprotein = ">Contig_1 1FFPMPQ*CCLFSSSYLPA*AIMVTGSCQTLVQDKV*CLASCHVMIPPCVFHPIYPLFTWSPLHPIISGCF*AVCMCSWIHALFVEIILSWSLSLPFTFEISSD"
protein = re.findall('M[A-Z]{6,}\*',myprotein)
print(protein)
['MVTGSCQTLVQDKV*', 'MIPPCVFHPIYPLFTWSPLHPIISGCF*']

```

---

To get the output you asked for while reading from a file, you can do:

```
import sys
import re

proteinFile = sys.argv[1]

with open(proteinFile, 'r') as f:
    for line in f:
        ## Split the long input line into the list 'sequences'
        sequences = line.rstrip('\n').split('>')

        for i, k in list(enumerate(sequences)):
            ## Skip the first element which will be empty.
            ## there must be a pythonic way of doing this, but I'm
            ## a Perl hacker, so I don't know it.
            if len(k) == 0:
                continue

            seqName, seq = k.split(' ')
            ## Get the frame
            frame = seq[0]
            ## Remove the frame from the sequence
            seq = seq[1:]
            # cleave sequence
            peptides = re.findall(r'(M[A-Z]{6,}\*)', seq)

            for idx, peptide in enumerate(peptides):
                print(">%s.%s.%d\n%s" % (seqName, frame, idx+1, peptide) )

```

---
Tags: fasta, python
---


# Question
Title: A Database for Plant Metabolome or Secondary Metabolites?
I am searching for a plant metabolome or secondary metabolite database for a college bioinformatics project. Up until now, I have found these two databases, as I am giving their links below;

and searched for different scientific articles and reviews;

But only one and updated database which could contain all known plants and all of their known metabolites would be perfect. Could you suggest me any database or any other source with a high coverage for my project?

# Answer
> 3 votes
PhytoMine is an InterMine interface to data from Phytozome, the Plant Comparative Genomics portal of the Department of Energy's Joint Genome Institute, which hosts 81 assembled and annotation genomes, from 78 viridiplantae species.

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


# Question
Title: FeaturePlot from Seurat: change its title
`FeaturePlot` is a function in Seurat package.

And in the vignette it is written that if we specify parameter `do.return = TRUE` it should return ggplot2 object. It is not working. My goal here is just to change the title of the plot. In case of violin plot I can do the following:

```
VlnPlot(object = seurat_object, features.plot = id, do.return = TRUE) + labs(title = endothelial_symbols[1])

```

But with `FeaturePlot` similar code fails to work:

```
FeaturePlot(object = seurat_object, features.plot = id, cols.use = c("grey", "blue"),
                          reduction.use = "tsne", do.return = TRUE) + labs(title = endothelial_symbols[1])

```

Giving the error:

> Error in FeaturePlot(object = seurat\_object, features.plot = id, cols.use = c("grey", : non-numeric argument to binary operator

Any suggestions would be greatly appreciated.

# Answer
> 3 votes
I tried with some data that I have and this is working for me:

```
p <- FeaturePlot(object = seurat_object, features.plot = id, cols.use = c("grey", "blue"),
                      reduction.use = "tsne", do.return = TRUE)

lapply(p, function(x){x + labs(title = endothelial_symbols[1])})

```

I think it is because FeaturePlot returns several ggplot objects according to the ids you put. This is returned as a list of plots so you have to iterate each one for labeling.

---
Tags: r, ggplot2, seurat
---


title: Hierarchical models with limma?

text:
# Question
Title: Hierarchical models with limma?
I have a dataset with (microarray) gene expression data that was sampled from the same individuals at multiple timepoints. Our exposure is a continuous variable, and because this was an observational study there is no consistent pattern of change in exposure over time (eg. one person's exposure may continuously increase, another's may go up and down, another may have no changes, etc.). The exposure is not expected to have any cumulative effect, so different timepoints within a person do not influence each-other. Our goal is to find genes where expression for the gene is associated with the level of exposure.

We are thinking of analyzing this as a hierarchical model, with timepoints nested within person (so a mixed model with a random effect for person). We will also need to control for a few covariates, some of which will be the same at all timepoints (race, sex) and some of which may change between timepoints (use of particular medications).

Is this something that it is possible to use limma for? I found that limma has a function called `duplicateCorrelation()` intended for block designs that some people have tried to use to fit mixed models, but I am not sure whether it can be used for this. Or if not limma, then what would you suggest?

# Answer
> 3 votes
Yes, you can use limma for this mixed model approach. Like you suggest, the random effect (persons) can be put in duplicateCorrelation().

Here is a similar example with RNAseq data, on bioconductor support site.

---
Tags: bioconductor, microarray, transcriptome, modelling
---


# Question
Title: differentially expressed intron analysis?
I have used Deseq for estimating the differentially expressed introns for the RNA seq data using the UCSC intron bed file. can Deseq be used for differentially expressed intron analysis?

# Answer
> 2 votes
Typical tools to use for this would be the DEXseq packages (as done here) or rMATS, which will calculate intron retention rates.

I would generally suggest not summing intronic counts by gene and using that in DESeq2 (I assume that's what you meant by DESeq, since there's no reason to use the original DESeq for anything). What you'd be getting if you did this is less expressed introns and more sample-specific bias in enriching for mature versus pre-mRNA.

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


# Question
Title: Behavior of `--reference` flag with samtools sort
I'm about to run samtools sort (version 1.6) and I'm reviewing the available configuration options.

```
$ samtools sort
Usage: samtools sort [options...] [in.bam]
Options:
  -l INT     Set compression level, from 0 (uncompressed) to 9 (best)
  -m INT     Set maximum memory per thread; suffix K/M/G recognized [768M]
  -n         Sort by read name
  -t TAG     Sort by value of TAG. Uses position as secondary index (or read name if -n is set)
  -o FILE    Write final output to FILE rather than standard output
  -T PREFIX  Write temporary files to PREFIX.nnnn.bam
      --input-fmt-option OPT[=VAL]
               Specify a single input file format option in the form
               of OPTION or OPTION=VALUE
  -O, --output-fmt FORMAT[,OPT[=VAL]]...
               Specify output format (SAM, BAM, CRAM)
      --output-fmt-option OPT[=VAL]
               Specify a single output file format option in the form
               of OPTION or OPTION=VALUE
      --reference FILE
               Reference sequence FASTA FILE [null]
  -@, --threads INT
               Number of additional threads to use [0]

```

I'm curious about the `--reference` flag. I've never used it before with `samtools sort` and I'm wondering what behavior it invokes. The online documentation doesn't even mention the flag. **Under what circumstances would I want/need to specify the reference genome?**

# Answer
You would only need to supply a reference when using CRAM output encoding. The CRAM format allows storing entries without explicitly encoding the base calls, since one can instead encode the differences to the reference genome.

With BAM output the `--reference` option has no effect.
> 4 votes

---
Tags: samtools
---


# Question
Title: tximport txidcol and geneidcol selecting wrong column
I'm trying to run a differential expression analysis of 2 samples of the same species, but with different treatments using DESeq2. I'm using Ensembldb to create a tx2gene dataframe for tximport, however tximport does not select the columns I specify, and instead always selects the tx\_biotype leading to the following output:

```
.                    25_quant    28_quant
nontranslating_CDS    10422.41    10022.74
protein_coding     22196722.59 24002269.26

```

I'm running tximport using the following params:

```
txi_tx <- tximport(files, type = "salmon", tx2gene = tx2gene,
                   txIdCol = "tx_id", geneIdCol = "gene_id")

```

I've gotten around this by deleting all columns from my tx2gene DF apart from tx\_id and gene\_id. However, I think this might be causing problems in the DESeq2 analysis somehow, and I would like to rule out that possibility.

This is what my `names(tx2gene)` looks like:

```
"tx_id"            "tx_biotype"       "tx_seq_start"     "tx_seq_end"
"tx_cds_seq_start" "tx_cds_seq_end"
"gene_id"          "tx_name"

```

What am I doing wrong?

# Answer
> 2 votes
By setting `type="salmon"` you're overriding your settings for things like `txIdCol`. The various `type` options (except `none`) are presets to set things correctly for the output from various programs like `salmon` or `kallisto`. If you need to specify custom column names, set `type="none"`.

As an aside, if I'm interpreting your question correctly you have only two samples in total. You can not perform a meaningful analysis with an unreplicated experiment, so if that is indeed your case then don't waste your time further with this dataset.

---
Tags: r, ensembl, tximport, tx2gene
---


# Question
Title: Off-by-one error when mapping Ensembl IDs to gene symbols
I have a matrix with `ensemblIds` which I am trying to convert to gene symbols using `BioMart`. Some of the `ensemblIds` will not be mapped, so I need to subset the rows of the matrix. It seems like `BioMart` adds one `ensemblId` or at least I see one more result ultimately than should be. Here is what I am trying to do:

```
library(dplyr)
library(tibble)
ids <- rownames(mtx_combined)
mart <- useDataset("mmusculus_gene_ensembl", useMart("ensembl"))
list <- getBM(filters= "ensembl_gene_id", attributes= c("ensembl_gene_id",
                           "mgi_symbol","description"),values=ids,mart= mart)
mtx_mapped <- as.matrix(as.data.frame(mtx_combined) %>%
rownames_to_column("gene") %>%
filter(rownames(mtx_combined) %in% list$ensembl_gene_id) %>%
column_to_rownames("gene"))

rownames(mtx_mapped) <- list$mgi_symbol   //Fails

```

Length of the returned `list` is `18337`, number of rows in the initial matrix is `18376`. After subsetting I expect to get `18337` rows, but getting one less - `18336`. What am I doing wrong here?

# Answer
There’s no 1:1 mapping between the IDs of different database schemas. As a consequence, some Ensembl gene IDs map to multiple MGI symbols, or to none (and vice versa).

Therefore, you can’t assume that the number of result rows will be equal to the number of query IDs.
> 6 votes

---
Tags: r, biomart
---


# Question
Title: Is there an efficient way to check an input BAM in R?
I'm writing a function in R for an R package which takes as input a BAM.

```
my_func = function(my_bam){
    if (!file.exists(my_bam)){
        stop("Input a valid BAM file.")
    }

    ## do awesome stuff with BAM

}

```

Trying to write robust software, I would like to do more of a check whether this BAM is actually a BAM. Above, I'm just using `file.exists()`.

I could also do something simple like "check whether the file has the extension `*.bam`" with

```
grepl('*.bam', my_bam)  ## TRUE for file.bam, FALSE for file.txt

```

However, is there a quick way to check that the input *is* actually a BAM? I would prefer to throw an error here, if the user inputs a \*txt or \*bedpe.

I don't think it's necessary to *validate* the BAM per se, whereby I believe most users rely upon Picard:

https://broadinstitute.github.io/picard/command-line-overview.html#ValidateSamFile

My first thought was Rsamtools for this operation:

https://bioconductor.org/packages/release/bioc/html/Rsamtools.html

However, it doesn't appear that there is a function here which checks the user inputs a BAM, e.g. `BamFile()` doesn't appear to throw an error

```
> library(Rsamtools)
> fake_bam = 'nonsense_string.bam'   ## a string I made up, not a real BAM
> BamFile(fake_bam)
class: BamFile
path: nonsense_string.bam
index: NA
isOpen: FALSE
yieldSize: NA
obeyQname: FALSE
asMates: FALSE
qnamePrefixEnd: NA
qnameSuffixStart: NA

```

How could I "check" that the user is inputting a BAM file? Perhaps pipe to samtools (which doesn't feel very stable in terms of software)?

# Answer
> 6 votes
You won’t know whether it’s a valid BAM file until you attempt reading it in full; even Devon’s method (which I recommend) only does a superficial check — i.e. it checks whether (parts of) the BAM *header* can be read correctly. However, this doesn’t guarantee that the file isn’t truncated or contains invalid data.

This is of course a more general problem, but it’s perhaps a bit stranger for BAMs (and other biological formats), as it’s relatively frequent to end up with technically invalid files that were produced by existing software. The only recourse is actually trying to use the file, and to handle errors appropriately.

For a *very quick check*™ you could also pass the file through `gunzip` and test whether the “magic bytes” at the beginning of the file equal `'BAM\1'`:

```
bgzf_stream = gzfile(filename, 'r')
on.exit(close(bgzf_stream))
magic = readChar(bgzf_stream, 4)
stopifnot(identical(magic, 'BAM\1'))

```

# Answer
> 3 votes
Take your `BamFile(fake_bam)` and pass it through `open()`. That will return a value if `fake_bam` is a real BAM (or CRAM or maybe SAM) file and through an error otherwise. You can then wrap all of that in a `trycatch()`:

```
trycatch(open(BamFile(fake_bam)), do something, do something else)

```

Alternatively, ditch the `trycatch()` and let `open()` return the error if an invalid file in entered.

# Answer
> 1 votes
Just expanding upon the answers by @DevonRyan and @KonradRudolph, here's one way to put this in an R package:

Either using `open(BamFile()` from Devon's answer, which needs `Rsamtools`:

```
my_func = function(my_bam){

    out = tryCatch(open(Rsamtools::BamFile(my_bam)), error = function(e) e)
    if (inherits(out, "error")){
        stop("Input a valid BAM file.")
    }

    ## do awesome stuff with BAM
}

```

Or from Konrad's answer

```
my_func = function(my_bam){

    magic = readChar(gzfile(filename, 'r'), 4)
    if (!identical(magic, 'BAM\1')){
        stop("Input a valid BAM file.")
    }

    ## do awesome stuff with BAM
}

```

---
Tags: r, sam, bioconductor, samtools
---


# Question
Title: 10X cellranger count, [error] The chemistry was unable to be automatically determined
While running `cellranger count` using the following `.slurm` file:

```
#SBATCH -p standard
#SBATCH -A overall
#SBATCH --time=48:00:00
#SBATCH --mail-type=end
#SBATCH --mail-user=nv4e@virginia.edu
#SBATCH --output=generate_count_%j.out
#SBATCH --error=generate_count_%j.err
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=20
#SBATCH --mem=120000

# Run driver program
FASTQS="/scratch/nv4e/kipnis/AHN3LLBGX5/outs/fastq_path/HN3LLBGX5"
TR="/scratch/nv4e/kipnis/refdata-cellranger-mm10-1.2.0"

cellranger count --id=17R_1 --transcriptome=${TR} --fastqs=${FASTQS} --sample=17R --expect-cells=2000

```

I got the error:

```
>2018-04-28 20:04:37 [runtime] (progress)
ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMI
STRY.fork0.chnk0: Indexing genome...
2018-04-28 20:04:49 [runtime] (progress)        ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMISTRY.fork0.chnk0: Building transcriptome...
2018-04-28 20:04:58 [runtime] (progress)        ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMISTRY.fork0.chnk0: Building kmer index...
2018-04-28 20:05:25 [runtime] (progress)        ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMISTRY.fork0.chnk0: Mapping reads...
2018-04-28 20:05:28 [runtime] (chunks_complete) ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMISTRY
2018-04-28 20:05:28 [runtime] (run:local)       ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMISTRY.fork0.join
2018-04-28 20:05:31 [runtime] (join_complete)   ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.CHEMISTRY_DETECTOR.DETECT_CHEMISTRY

2018-04-28 20:05:34 [runtime] (ready)           ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.SETUP_CHUNKS
2018-04-28 20:05:37 [runtime] (split_complete)  ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.SETUP_CHUNKS
2018-04-28 20:05:37 [runtime] (run:local)       ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.SETUP_CHUNKS.fork0.chnk0.main
2018-04-28 20:05:40 [runtime] (failed)
ID.17R_1.SC_RNA_COUNTER_CS.SC_RNA_COUNTER.SETUP_CHUNKS

>[error] The chemistry was unable to be automatically determined. This can happen if not enough reads originate from the given reference. Please verify your choice of reference or explicitly specify the chemistry via the --chemistry argument

```

What is going wrong here? How could I correct it?

# Answer
The issue was that I supplied wrong indexes during `fastq` files generation, however, they were generated successfully, and only the next step failed with not that informative error.
> 2 votes

---
Tags: 10x-genomics, cellranger, slurm
---


# Question
Title: How to perform functional analysis on a gene list in R?
From an RNA-seq experiment I have about 17000 gene ids for 2 sample conditions arranged according to their log2 fold changes when compared to a control. I need to annotate these, but I've never done annotation before and am wondering how to do this in R.

I'm interested in up and down regulation of the pathways of the involved genes. I'd love to get enrichment and p-value analysis too. I'm interested in whether or not a package similar to PANTHER exists in R for this analysis. I've been reading how different annotation programs can be better or worse; it seems PANTHER is much better than DAVID, so I was wondering if there is a package that provides PANTHER-like analyses in R.

There seems to be multiple packages available, are there any that stand out as being the best?

I'm primarily interested in human samples and annotated pathways.

# Answer
> 7 votes
N.B., I'm avoiding discussion of "best", since that's more or less impossible to answer.

Your question can actually be divided into two:

1. What's a good tool for pathways analysis (ideally in R)?
2. What are good sources of pathway information?

For (1), there are a few different possibilities, but I prefer either roast or camera from limma or goseq (a stand-alone bioconductor package). Of those, I expect what roast or goseq are doing are the most similar to what panther is doing in the sense that they're not competitive. However, since genes do correlate with each other, I prefer competitive tests like that offered by camera.

For (2), your options are basically KEGG or reactome. KEGG is problematic, since you need a license for anything remotely recent (there are some R packages to get around this, but I haven't a clue how kosher they are legally). Given that, reactome might be your best bet.

Having written all of that, the best database out there is what IPA uses. This is a commercial product (I have no affiliation), but you can get a demo license that will work for a couple weeks at least. If you're doing analyses relatively often then it makes sense to spend a bit of cash on a license (maybe share with neighboring groups).

Regardless of what you use, you might also want to use pathview, which can produce easily interpretable plots related to pathway enrichment.

# Answer
> 4 votes
You requested a tool similar to PANTHER but in R. First, the PANTHER (Protein ANalysis THrough Evolutionary Relationships) tool does a classification based on evolutionarily related proteins, gene ontologies ( molecular function, and biological process) and pathways.

AFAIK there isn't a tool in R that integrates all these into one, but several packages do pieces of it. However, my observations after using it a couple of times is that is based a lot in gene ontologies. Gene Ontologies (GO) are not pathways, but the best tool I found for GO are topGO (A bit hard to use but useful to get insights on the biology) which test over representations of a term in a given list.

# Answer
> 1 votes
I would go for the topGO package. The logical steps in your pipeline would be:

1) Build a matrix with the following variables:

|gene\_ID|GO\_ID|GO\_term\_description| You can enter a list of uniprot IDs to uniprotkb and retrieve the GO annotations for each gene in your list/file

2) filter that list for the genes that were captured in your differential expression analysis even tho they are not differentially expressed. This filter is defined by DESeq/Sleuth by "at least 4 read counts in at least 49% of the libraries"

3)Define this gene list as gene\_universe

4) define your gene set and get the gene\_IDs and if possible, also a ranked score (p-val/FDR) of the genes you are interested in

5) plot the graphs and get tables of enriched GO terms (these would be ranked by p-value without need of performing correction analysis

Note: you are going to need some data manipulation functions; I would read the "tidyverse" documentation

---
Tags: r, sequence-annotation
---


title: Subset on multiple genes in Seurat

text:
# Question
Title: Subset on multiple genes in Seurat
I know that I can do subsetting on just one gene in `Seurat`:

```
seurat_subset <- SubsetData(seurat_object, subset.name = neuron_ids[1], accept.low = 0.1)

```

However, I want to subset on multiple genes. Is there a way to do that? I just do not want to do manual subsetting on 10 genes, then manually getting `@data` matrix from each subset, and recreating seurat object afterwards.

# Answer
I was able to achieve this in the following way:

```
require(data.table)
cell_names <- vector(mode="character")
for (i in 1:length(neuron_ids)) {
  seurat_subset <- SubsetData(seurat_object, subset.name = neuron_ids[i], accept.low = 0.1)
  genes <- colnames(seurat_subset@data)
  cell_names <- c(cell_names, genes[! genes %chin% cell_names])
}
seurat_subset <- SubsetData(seurat_object, cells.use = cell_names)

```

Would be interesting to know if `Seurat` provides such functionality out of the box.
> 6 votes

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


# Question
Title: Does the galGal5 chicken assembly have a chromosome 29?
The chromosome sizes at UCSC don't seem to contain chr29:

ftp://hgdownload.soe.ucsc.edu/goldenPath/galGal5/bigZips/galGal5.chrom.sizes

It has a chr28 and a chr30. Am I missing something or is there some piece of history that led to the omission of this chromosome?

# Answer
The current ensembl entry doesn't have a 29 either. The archived ensembl assembly lacks 29 30, and 31 and 33 and LGE64.

The chromosomes after 30 are tiny, so they might not be visible in a karyotype. They probably realized that "chr 29" was really attached to some other chromosome.
> 2 votes

---
Tags: genome, assembly, reference-genome, ucsc
---


title: BioMart Ensembl id imperfect retrieval

text:
# Question
Title: BioMart Ensembl id imperfect retrieval
I have a list of gene symbols:

```
c("cd45", "Tmem119", "CD11b", "P2Yr12", "Csf1r", "Bst2", "Cd74",
"Cx3cr1", "Trem2", "Lyz2", "GLAST", "GFAP", "ALDH1L1", "Aldoc", "Aqp4",
"Glul", "S100a", "Olig1", "Olig2", "Olig3", "Mbp", "Pdgfra", "Pecam",
"Cldn5", "Cldn10", "Epas1", "Crip1")

```

If I feed it to `BioMart` to get `ensembl ids`:

```
mart <- useDataset("mmusculus_gene_ensembl", useMart("ensembl"))
list <- getBM(filters= "mgi_symbol", attributes= c("ensembl_gene_id",
                           "mgi_symbol","description"),values=symbols,mart= mart)

```

the following genes are missed:

```
"cd45"   "cd11b"  "p2yr12" "glast"  "s100a"  "pecam"

```

All of them are pretty well-known genes, and I can manually find their `ensembl ids` by googling it, for instance:

http://www.informatics.jax.org/marker/MGI:97810

I tried supplying aliases, but it does not change the output. So, from my understanding `BioMart` is not working properly either because I am doing something wrong here or because `BioMart` itself is not a good tool to use. Is there a better way of getting the mapping that would map all of the gene symbols?

# Answer
> 4 votes
The genes that are missed are probably not official mgi symbols.

You might wanna look them up at mgi:

```
cd45 -> Ptprc

cd11b -> Itgam

p2ry12 -> P2ry12 (?) (CAPS sensitive?)

glast -> Slc1a3

s100a -> S100a1

pecam -> Pecam1 (?)

```

I suspect biomaRt is not the problem here, the only better way I can come up with is to download the ENSEMBL annotation (GTF) file, and get the symbol out of this annotation with AWK or in R.

Here an example in R, since you don't want to use AWK:

```
file <- read.table("Mus_musculus.GRCm38.78.gtf", sep = "\t", comment.char = "#")

exons <- subset(file, V3 == 'exon')

V9 <- file[,"V9"]

ensembl_names <- sub("gene_id ","",V9)

split <- strsplit(ensembl_names, ";")

splitting <- sapply(split, function(x) x[1])

split_symbol <- strsplit(ensembl_names, "gene_name ")

splitting_symbol <- sapply(split_symbol, function(x) x[2])

splitting_symbol2 <- strsplit(splitting_symbol, "; ")

symbol_names <- sapply(splitting_symbol2, function(x) x[1])

ensemble_symbol_big <- cbind(splitting, symbol_names)

ensembl_symbol <- unique(ensemble_symbol_big)

colnames(ensembl_symbol) <- c("ENSEMBL", "Symbol")

write.table(ensembl_symbol, "ensembl_symbol.txt", sep = "\t", row.names = F)

```

---
Tags: r, biomart
---


title: Average of logCPM

text:
# Question
Title: Average of logCPM
I am a bit in doubt how to calculate correctly the mean (average) in the logarithmic scale. I know e.g., that a division in the normal scale becomes a subtraction in the logarithmic scale, and so forth, that's why I was wondering what the correct way is to calculate the mean (since the mean is calculated by additions and divisions in the normal scale).

So for example, if I have in RNAseq data a group with 3 samples, and have log2 CPM values for each sample. How can I calculate the mean of the group correctly?

Do I first need to get CPM (without log2), calculate the mean and then apply log2? Or can I just ignore the fact that the log2 CPM values are in the log scale, and just calculate the mean as if it was the normal scale (add all up and divide by *n*)?

# Answer
> 2 votes
The mean of log-transformed values is the log-transformed geometric mean of the untransformed values (i.e., $log\_{2}(geometric~mean)$). This can be useful when the values you want to summarize have vastly different ranges (so one doesn't dominate the average). I would expect that to generally **NOT** be the case for CPMs of the same gene, so it's unclear to me what the point of that would be. It can also be difficult to interpret these results due to the geometric mean being dependent upon the variance:

```
> originalValues = rnorm(1000, mean=10, sd=c(rep(1, 500), rep(2, 500)))
> mean(originalValues[1:500])
[1] 9.991464
> mean(originalValues[501:1000])
[1] 10.04407
> log2(mean(originalValues[1:500]))
[1] 3.320696
> log2(mean(originalValues[501:1000]))
[1] 3.328272
> mean(log2(originalValues[1:500]))
[1] 3.313263
> mean(log2(originalValues[501:1000]))
[1] 3.300103

```

Above I created two groups of 500 values each. Each group has approximately the same mean (10) but with different variance (1 vs 2). As you can see, the second group of 500 values randomly has a slightly higher average. The log2(mean) preserves this difference, but the log2(geometric mean) inverts the relationship such that group 2 now appears to have a lower mean value. This isn't really wrong, it's the difference between mean and geometric mean, but it's something that needs to be kept in mind.

This has implications if your want to compare values between genes, since your comparison metric is then not just dependent upon the mean, but also the variability of the data (such that more variable genes will have systematically lower geometric means).

In short, the geometric mean is useful when values are on different scales (such as CPM values of all genes in a sample, but probably not CPM values of the same gene between samples). If that's not the case, then you're mostly removing interpretability from the results.

---
Tags: rna-seq
---


# Question
Title: How to get rid of the temp files while using "featureCounts" for extracting readcounts from bam files?
I used this command for extracting read counts for almost 60 samples using featureCounts.

```
featureCounts -a /gencode.v27.primary_assembly.annotation_nochr.gtf -F GTF -p -s 2 -T 8 -o counts.txt sample.sorted.bam

```

I see many "temp-core-100373-A4BF0116E148.sam-TH00-BK004414.tmp" files like this dumped in my directory. I'm able to delete them only manually. Is there a way to get rid of these tmp files from my directory while using the command?

# Answer
You shouldn't delete them while featureCounts is running, it's creating read-name-sorted temporary files so it can more easily look at both mates together. Just wait for it to finish running and it will clean up the files itself. Doing so beforehand is liable to cause an error.
> 2 votes

---
Tags: rna-seq, sam, rna, reads, gtf
---


title: get gene lines from gtf file

text:
# Question
Title: get gene lines from gtf file
I would like to retrieve gene lines from a GTF file for which I only have exons & transcripts lines (output from Cufflinks) and alternative splicing possible. I need gene lines for compatibility with a pipeline dealing with GTF in Ensembl format.

Ideally, I would have gene lines representing the longest transcript - but I am open to discussion about that. The aim would be to separate the genome in genic & intergenic portions - boundaries being defined by start-end coordinates of the genes.

Sample input (e.g. gene `"CUFF.105"` with 5 transcripts):

```
cat subset_105.test | awk '$3=="transcript"{print $0}'

tig00000046 Cufflinks   transcript  26170   42766   202 +   .   gene_id "CUFF.105"; transcript_id "CUFF.105.3"; FPKM "0.2320304094"; frac "0.197788"; conf_lo "0.164885"; conf_hi "0.299176"; cov "5.188855";
tig00000046 Cufflinks   transcript  26170   42766   266 +   .   gene_id "CUFF.105"; transcript_id "CUFF.105.1"; FPKM "0.3061260755"; frac "0.249779"; conf_lo "0.228860"; conf_hi "0.383392"; cov "6.845843";
tig00000046 Cufflinks   transcript  26170   39469   470 +   .   gene_id "CUFF.105"; transcript_id "CUFF.105.2"; FPKM "0.5403628578"; frac "0.161731"; conf_lo "0.372512"; conf_hi "0.708214"; cov "12.084038";
tig00000046 Cufflinks   transcript  28928   39469   1000    +   .   gene_id "CUFF.105"; transcript_id "CUFF.105.4"; FPKM "1.1485378983"; frac "0.233911"; conf_lo "0.853387"; conf_hi "1.443688"; cov "25.684547";
tig00000046 Cufflinks   transcript  29614   42766   181 +   .   gene_id "CUFF.105"; transcript_id "CUFF.105.5"; FPKM "0.2087076495"; frac "0.156792"; conf_lo "0.137112"; conf_hi "0.280304"; cov "4.667292";

```

What it came to my mind is:

* sort the file by `$4` (numeric: n) and `$5` (reverse: nr) (`sort -k4,4n -k5,5nr`)
* take the 1st line
* remove unnecessary fields from `$9`

But this would fail for transcripts annotated in the reverse strand, and I won't be able to deal (how to fill) the `$6` (score - a floating point value).

Any ideas of an existing method that could help me retrieve those gene lines?

Options to retrieve those lines directly from Cufflinks would also be accepted (I didn't find a suitable parameter in the manual).

# Answer
Assuming none of the genes have duplicate entries across chromosomes (or strand or other biologically implausible things) in R one can:

```
library(rtracklayer)
library(GenomicRanges)
gtf = import.gff("file.gtf")
grl = split(gtf, gtf$gene_id)  # This produces a GRangesList
grl = endoapply(grl, function(x) {
    foo = x[1]
    foo$type = 'gene'
    start(foo) = min(start(x))
    end(foo) = max(end(x))
    return(c(foo, x))
    })
gr = unlist(gr)

```

You can then either write `gr` out to a file (either with `export()` or just printing the rows).
> 1 votes

---
Tags: sequence-annotation, file-formats, gtf
---


# Question
Title: Which tools for differential expression analysis in scRNA-Seq?
I am starting to run analysis for differential expression in scRNA-Seq. Which tools are available for this kind of analysis?

Can tools for bulk RNA-Seq like DESeq be used for scRNA-Seq?

# Answer
According to Jaakkola et al. (ref. 1):

> Based on our comparisons, DESeq and Limma without any modifications are not suitable for scRNA-seq data analysis, and yet, they have performed well in the context of bulk RNA-seq data.

They also provided a summary of the comparison between 5 tools (2 designed for Single-Cell, 2 designed for bulk, 1 general):

In my case I ended up using MAST (ref. 2, github).

1. Jaakkola, M. K., Seyednasrollah, F., Mehmood, A. & Elo, L. L. Comparison of methods to detect differentially expressed genes between single-cell populations. Briefings in Bioinformatics 18, 735–743 (2017).
2. Finak, G. et al. MAST: A flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data. Genome Biol. 16, 1183 (2015).
> 1 votes

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


# Question
Title: How to validate that BAMs have been downloaded correctly?
I currently have several hundred BAM files which were downloaded by someone else. These have remained untouched---before working with them, I would like to double-check that these BAMs have been fully downloaded.

I don't MD5 Checksums to look at.

In order to validate the BAMs, I would use something like Picard's ValidateSamFile, https://broadinstitute.github.io/picard/command-line-overview.html#ValidateSamFile

```
java -jar picard.jar ValidateSamFile \
      I=input.bam \
      MODE=SUMMARY

```

I could write a bash script, and run through these BAMs, and then read through the Picard output. That is checking whether the BAM is valid though, not whether it was downloaded fully (e.g. the BAM could be downloaded correctly, but be inherently flawed in some other way).

If the BAM is downloaded fully, isn't there some "end of file" bytes I could quickly check? Are there other more efficient methods?

# Answer
> 9 votes
`samtools quickcheck` is all you need. From the manual:

> Quickly check that input files appear to be intact. Checks that beginning of the file contains a valid header (all formats) containing at least one target sequence and then seeks to the end of the file and checks that an end-of-file (EOF) is present and intact (BAM only).
>
> Data in the middle of the file is not read since that would be much more time consuming, so please note that this command will not detect internal corruption, but is useful for testing that files are not truncated before performing more intensive tasks on them.
>
> This command will exit with a non-zero exit code if any input files don't have a valid header or are missing an EOF block. Otherwise it will exit successfully (with a zero exit code).

---
Tags: sam, data-download, picard, validation
---


title: How to enlarge a section of y-axis?

text:
# Question
Title: How to enlarge a section of y-axis?
I am using `ggplot2` to plot a scatter plot.

```
library(ggforce) # required by facet_zoom

g <- ggplot(data.plot, aes(x = Methyl_Average, y = Average_Gene_Expression))
g1 <-  g + geom_point(color = "steelblue") + theme_bw(base_family = "Times") +
       ggtitle( "Scatter Plot") +
       theme_bw(base_family = "Times") +
       theme(plot.title = element_text(hjust = 0.5))

```

As most of my data points are in the range of 0 - 200 on y-axis, I want to zoom this section. I tried doing it using `facet_zoom` from `ggforce` package :

```
g1 + facet_zoom(y = Average_Gene_Expression > 0 & Average_Gene_Expression < 200)

```

Instead of two plots in one figure, I want one plot (like first figure), where 80% of y-axis displays the data points in the range of 0 - 200 and 20% of the y-axis displays the remaining and x-axis unaltered. Which function/package should I use do it in R?

# Answer
I would strongly discourage you from making discontinuous axis, it's going to be very confusing for a reader.

The facet plot you proposed seems like a good solution to me. Alternatively you can use log transformation. To demonstrate I made it on simulated data that look appox like yours :

```
set.seed(940401)

data.plot <- data.frame(Methyl_Average = c(0.05, rbeta(699, 2, 5) + 0.15, runif(300, 0.18, 0.85)),
                        Average_Gene_Expression = c(165, exp(rnorm(999, 3))))

library(ggplot2)
require(gridExtra)

g <- ggplot(data.plot, aes(x = Methyl_Average, y = Average_Gene_Expression))
g1 <-  g + geom_point(color = "steelblue") + theme_bw(base_family = "Times") +
       ggtitle( "Scatter Plot") +
       theme_bw(base_family = "Times") +
       theme(plot.title = element_text(hjust = 0.5))

sc_plot <- g1 +
           labs(x = "Differential Methylation", y = "Differential Expression")
logsc_plot <- g1 +
              labs(x = "Differential Methylation", y = "Differential Expression (log10 scale)") +
              scale_y_continuous(trans='log10')

grid.arrange(sc_plot, logsc_plot, ncol = 2)

```
> 5 votes

---
Tags: r, visualization, ggplot2
---


# Question
Title: Demultiplexing and preprocessing for custom single-nucleus Drop-seq data
I am trying to reproduce the preprocessing of paired-end sequencing reads in Lake et al. (ref. 1).

> First, paired-end reads were filtered out if read 1 had more than four non-T bases in the last ten bases (to remove all non-poly(T)-captured contaminated reads), or had one or more bases with a poor quality score (\<10). Cell barcode and UMI information were then inferred from the first 12 bases and the next 8 bases of read 1, respectively. The right mate of each read pair was trimmed to remove any portion of the SMART adaptor sequence or large stretches of poly(A) tails (6 consecutive bp or larger).

I was able to filter out reads as described in the first sentence.

However, I am unsure on which tool or procedure to use to demultiplex based on cell barcode and UMI, and to trim the right mate. Do you know what I can use or something I can read to understand how to do it?

1. Lake, B. B. et al. Integrative single-cell analysis of transcriptional and epigenetic states in the human adult brain. Nature Biotechnology 36, 70–80 (2017).

# Answer
I didn't notice the sentence:

> Paired-end sequencing reads were processed largely as described (http://mccarrolllab.com/wp-content/uploads/2016/03/Drop-seqAlignmentCookbookv1.2Jan2016.pdf ), with additional correction steps.

This should provide instructions on how to lead the preprocessing and mapping of data.
> 0 votes

---
Tags: scrnaseq, reads, trimming, demultiplexing
---


# Question
Title: Transform and feed data into DESeq2 with DESeqDataSetFromMatrix
There is a `normalized` expression matrix. I split it into two and want to do `DE` on the two cells' subsets. I am having trouble transforming it into the format that `DESeq2` would accept. The function that I would think I need to use is the following:

```
dds <- DESeqDataSetFromMatrix(countData = cts,
                          colData = coldata,
                          design= ~ batch + condition)

```

It would be perfect if I could somehow feed two matrixes separately, but I can not find it. Any suggestions would be greatly appreciated.

> Update

I tried the following:

```
cells <- c(cell_names, cell_names_cyt)
groups <- c(rep("NE", times=length(cell_names)), rep("EX",
times=length(cell_names_cyt)))

colData <- data.frame(
  Cells <- cells,
  Groups <- groups
)

dds <- DESeqDataSetFromMatrix(countData = df,
                          colData = colData,
                          design= ~ Groups)

```

But getting the error:

> Error in DESeqDataSet(se, design = design, ignoreRank) : some values in assay are not integers

Of course, they are floats since I already normalized the matrices.

> Update

I did the following since I am using `Seurat`:

```
df <- as.data.frame(as.matrix(seurat_object@raw.data))
df <- df[,colnames(df) %in% cells]

```

But then the following:

```
dds <- DESeqDataSetFromMatrix(countData = df,
                          colData = colData,
                          design= ~ Groups)

```

Gave me an error:

> Error in DESeqDataSet(se, design = design, ignoreRank) : all variables in design formula must be columns in colData

# Answer
> Error in DESeqDataSet(se, design = design, ignoreRank) : all variables in design formula must be columns in colData

I think that means just what it says...you don't have a "Groups" column in coldata. Did you spotcheck your ColData data frame to make sure it looks the way it should?

```
> cells <- c("1","2","3")
> groups <- c("a", "b", "c")
> test <- data.frame(Cells <- cells, Groups <- groups)
> test
  Cells....cells Groups....groups
1              1                a
2              2                b
3              3                c
> colnames(test)
[1] "Cells....cells"   "Groups....groups"

```
> 1 votes

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


title: update dbSNP ID

text:
# Question
Title: update dbSNP ID
I have used dbSNP build 138 to tag 'rs' ID in my VCF file. But later I came to know that now dbSNP build 150 has been released. I want to update the rs IDs in my VCF file.

I tried but I found a problem for which I need suggestion. With dbSNP build 138, I am getting 450076 total variants, 248236 as known variants and the rest as novel. But when I use the latest dbSNP file, `All_20170710.vcf.gz`, I am getting same total variants but 0 as known and all as novel which is definitely wrong.

When I checked the reason for this, I found that in the `All_20170710.vcf.gz` VCF file, all chromosomes are written as numbers like `1`, `2`, instead of `chr1`, `chr2` while in my vcf file, all chromosomes are written as `chr1`, `chr2`. The dbSNP build 138 file also has the chromosome names with the `chr` prefix, and that's why 138th build was giving good results.

One solution would be add the `chr` word before chromosome number. Is that an efficient way? How do I do that? If not, what is the correct way to use the latest dbSNP `All_20170710.vcf.gz` to tag rs IDs in a vcf file which contains ~32 Billion SNPs instead of `dbsnp_138.vcf` which contains only ~6 Billion SNPs? Any suggestions??

Thanks

# Answer
> 2 votes
It sounds like you downloaded the dbSNP release for the GRCh38 or GRCh37 genome releases from the Genome Reference Consortium (GRC) which name the chromosomes without the `chr` prefix and you were previously using a file that had been built using either hg19 or hg38 which use the `chr`.

As far as I know, dbSNP only offers vcf files built against the GRC genome builds, so you (or someone else) had probably already modified the old file you were using and changed the chromosome names.

Also note that if you want the latest version, you should download ftp://ftp.ncbi.nih.gov/snp/organisms/human\_9606/VCF/All\_20180418.vcf.gz.

Anyway, yes, you can simply change the VCF to use `chrN` instead of `N`. Note, however, that you will also need to change the name of the mitochondrial chromosome from `MT` to `chrM`

You can do all this very easily with `sed`:

```
zcat All_20170710.vcf.gz |
    sed -E '/^[^#]/{s/^([^#]+)/chr\1/; s/^chrMT/chrM/}' > All_20170710.fixed.vcf

```

The command above will work on most `sed` implementations, including Linux (GNU sed) and also macOS (BSD sed).

---
Tags: vcf, sequence-annotation, snp, ngs
---


# Question
Title: Which are the use cases for the methods for DE in Seurat
In `Seurat` we can specify multiple methods for finding `DE` genes.

I am wondering when should we use which one, the use cases. As I understand, whole `Seurat` package is used for single cell sequencing analysis. I ran different methods for my dataset that was produced from brain cells, around `800` cells in each sample, by single cell sequencing. Now I am confused which methods are the most appropriate to use, which results are the most reliable. For single cell sequencing when would you prefer to use one method over the other?

# Answer
You can take a look at the recently published article: Bias, robustness and scalability in single-cell differential expression analysis.

> We evaluated 36 approaches using experimental and synthetic data and found considerable differences in the number and characteristics of the genes that are called differentially expressed. Prefiltering of lowly expressed genes has important effects, particularly for some of the methods developed for bulk RNA-seq data analysis. However, we found that **bulk RNA-seq analysis methods do not generally perform worse** than those developed specifically for scRNA-seq.

Then you can decide which method fits better to your type of scRNA-seq protocol. Since Seurat offers DESeq2 could be a good idea to use it, but this depends a lot on your data.
> 5 votes

---
Tags: r, single-cell, seurat
---


# Question
Title: A reliable fetcher of short read using SRA/ENA accession
I am trying to build a workflow that gets data automatically from databases of sequencing reads (SRA - sequence read archive, ENA - european nucleotide archive). Till now I was pulling everything from SRA and my script have worked, but now I encountered a first entry that is only in ENA (ERR2135453), not in SRA. And I found that my download script has its limitations.

My bash script looks like this (requires aspera `ascp` and `fastq-dump`) and I wanted to make changes that are commented now

```
#!/bin/bash

ACCESION=$1
OPENSSH=/my/part/to/key/asperaweb_id_dsa.openssh

# if [ ${ACCESION:0:3} == SRA ]; then
    URL=anonftp@ftp-private.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/${ACCESION::3}/${ACCESION::6}/"$ACCESION"/"$ACCESION".sra data/
# else
#    URL=ena_url
# fi

ascp -QT -i $OPENSSH $URL

fastq-dump data/"$ACCESION".sra --outdir data/ --split-files --gzip
rm data/"$ACCESION".sra

```

The problem is in getting ENA url. Unlike SRA url that I can easily get from the accession (`anonftp@ftp-private.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/${ACCESION::3}/${ACCESION::6}/"$ACCESION"/"$ACCESION".sra data/`), the ftp server of ENA is bit more confusing. The url corresponding to ERR2135453 accession is `ftp://ftp.sra.ebi.ac.uk//vol1/err/ERR213/003/ERR2135453`. Where `003` there comes from? Also not all reads are in `err` subdirectory.

Is there an easy way how to download reads from ENA knowing just the accession number?

--- edit ---

I found on the relevant ENA site about downloading data using aspera that ftp urls should be possible to derive only from accession also in ENA.

> Submitted read data files are organised by submission accession number under vol1/ directory in ftp.sra.ebi.ac.uk:
>
> ftp://ftp.sra.ebi.ac.uk/vol1//
>
> where contains the first 6 letters and numbers of the SRA Submission accession.
>
> For example, the files submitted in the SRA Submission ERA007448 are available at: ftp://ftp.sra.ebi.ac.uk/vol1/ERA007/ERA007448/.

The reason why there is this strange `003` will probably be that it's quite recent submission and it has not been properly categorized yet (other accessions I tried were on "expected" urls).

However, the problem still persist. Is there a way how to fetch data of new submissions using accessions only?

# Answer
> 1 votes
Apparently there is no **100% reliable** solution, but in the most cases, you should be able to use the enaBrowserTools python scripts (enaDataGet) to quite reliably download data.

The run ERR2135453 is a special case because the sample submission (ERS1939590) for it has been cancelled. The solution for these reads is to download the data manually from the ftp location.

This answer is based on kind assistance of Suran from ebi helpdesk.

# Answer
> 1 votes
Here are two python functions I use for getting stuff from ENA based on an accession:

```
import os
import itertools
from urllib.request import urlopen

def fetch_ENA_files(accession):
    """Get the names of the files matching the ENA accession"""

    # Get the paths to the files
    url = "http://www.ebi.ac.uk/ena/data/warehouse/filereport?accession=%s&result=read_run&fields=fastq_ftp" \
          % accession

    try:
        paths = urlopen(url).readlines()[1:]
    except:
        E.debug("couldn't access %s" % url)
        raise

    paths = list(itertools.chain.from_iterable(
        [p.strip().split(";") for p in paths]))
    filenames = [os.path.basename(x) for x in paths]
    dl_paths = [os.path.join(*x.split("/")[1:]) for x in paths]
    return filenames, dl_paths

def fetch_ENA(dl_path, outdir, protocol="ascp"):
    """Fetch fastq from ENA given accession"""

    if protocol == "ascp":
        statement = """ascp -QT -P33001 -l 100M -i $ASCP_KEY_PATH
                       era-fasp@fasp.sra.ebi.ac.uk:/%(dl_path)s %(outdir)s """ % locals()
    elif protocol == "http":
        fn = os.path.basename(dl_path)
        outfile = os.path.join(outdir, fn)
        statement = "wget -O %(outfile)s ftp://ftp.sra.ebi.ac.uk:/%(dl_path)s" % locals()

    return statement

```

`fetch_ENA_files` queries the ENA database for the names of all the files associated with a accession and returns a tuple of lists. The first entry in the tuple contains the filenames (which you might need later in your application), the second contains a list of paths on the ENA servers where those files can be found.

`fetch_ENA(dl_path, outdir, protocol)` takes a download path and downloads it to `outdir` using with ftp or ascp.

You might use it like:

```
filenames, paths = fetch_ENA_files("ERR2135453")
for dlp in paths:
    bash_statement = fetch_ENA(dlp, "downloads", "ascp")
    os.system(bash_statement)

```

(we don't we actaully use our modules to send it to the cluster)

---
Tags: data-download, sratoolkit
---


# Question
Title: GATK CombineVariants complains the contig order in the VCF files
I have called variants on two strains of *C. elegans* separately. I now want to merge the VCF files into one using the following code:

* Create a sequence dictionary of the reference sequence
* Sort the VCF files with Picard
* Merge the sorted VCF files using GATK

```
picard CreateSequenceDictionary \
    REFERENCE=c_elegans.PRJNA13758.WS263.genomic.fa \
    OUTPUT=c_elegans.PRJNA13758.WS263.genomic.dict

picard SortVcf INPUT=strain1.vcf \
    OUTPUT=strain1sorted.vcf \
    SEQUENCE_DICTIONARY=c_elegans.PRJNA13758.WS263.genomic.dict

picard SortVcf INPUT=strain2.vcf \
    OUTPUT=strain2sorted.vcf  \
    SEQUENCE_DICTIONARY=c_elegans.PRJNA13758.WS263.genomic.dict

GATK --analysis_type CombineVariants \
    -R c_elegans.PRJNA13758.WS263.genomic.fa \
    --variant strain1sorted.vcf \
    --variant strain2sorted.vcf \
    -o all.vcf \
    -genotypeMergeOptions UNIQUIFY

```

The last command gives me the following error message:

```
ERROR MESSAGE: Input files variant and reference have incompatible contigs. Please see https://www.broadinstitute.org/gatk/guide/article?id=63 for more information. Error details: The contig order in variant and reference is not the same; to fix this please see: (https://www.broadinstitute.org/gatk/guide/article?id=1328),  which describes reordering contigs in BAM and VCF files..
##### ERROR   variant contigs = [I, II, III, IV, MtDNA, V, X]
##### ERROR   reference contigs = [I, II, III, IV, V, X, MtDNA]

```

But I have sorted the VCF files using Picard, so I don't know what else to do.

Your help is appreciated.

Update: I asked this question on the GATK forum on 4th May 2018 (today).

# Answer
> 2 votes
GATK cannot combine VCF files generated by Sniffles, the variant caller that I used to call structural variants (answer from the GATK team).

# Answer
> 3 votes
Answered in the GATK forum

> Try regenerating the index files for your VCFs. Picard SortVcf doesn't do it for you iirc, so when GATK looks at the index files (before opening the VCFs themselves) it sees the old order and complains.

---
Tags: vcf, variant-calling, gatk, picard, variation
---


# Question
Title: Generating multiple NCBI URLs with eutils
I would like to create a URL that takes input from different arrays to generate a search. The data I would like to retrieve is from NCBI. Below you can find the code I wrote. But I can't get it to work; I am stuck. Thank you in advance.

I should have something like this as result but i got nothing. https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=gene&retmode=xml&retmax=10&sort=relevance&term=A20%20AND%20Homo%20sapiens%20%5borgn%5d%20AND%20alive%5bprop%5d

```
$eutils = "http://eutils.ncbi.nlm.nih.gov/entrez/eutils";
$db = "gene";
$gene_symbols = array("A20","HBB","ABO");
$species = array("homo sapiens", "Escherichia coli","Saccharomyces cerevisiae");
$links = array();
foreach ($gene_symbols as $key => $GS){
    $SP = $species[$key];
    $query = "$GS AND $SP [orgn] AND alive";
    $esearch = "$eutils/esearch.fcgi?db=$db&retmode=xml&retmax=10&sort=relevance&term=";
    #echo $esearch;
    $links[] = get($esearch.$query);
}
foreach($links as $lk){
    xml = simplexml_load_file($lk) or die("feed not loading");
}
var_dump(xml);

```

As output I got nothing.

But I would like to have this: https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=gene&retmode=xml&retmax=10&sort=relevance&term=A20%20AND%20Homo%20sapiens%20%5borgn%5d%20AND%20alive%5bprop%5d And this for every gene. For some reason it put them all after each other and the foreach loop doesn't seem to work. I got no errors.

I run this code in php using visual studio code. The result should be visualised in the webbrowser (http://localhost/test.php).

# Answer
I don't know PHP, so I'm sure there are better ways of doing this, but you had various issues in your code. This seems to be a working version of what you're after:

```
<?php
$eutils = "http://eutils.ncbi.nlm.nih.gov/entrez/eutils"; $db = "gene"; $gene_symbols = array("A20","HBB","ABO");

$species = array("homo%20sapiens", "Escherichia%20coli","Saccharomyces%20cerevisiae");
$links = array();

foreach ($gene_symbols as $key => $GS){
  $SP = $species[$key];
  $query = "$GS%20AND%20$SP%20[orgn]%20AND%20alive";

  $esearch = "$eutils/esearch.fcgi?db=$db&retmode=xml&retmax=10&sort=relevance&term=";

  echo "$esearch$query\n";
  $links[] = file_get_contents($esearch.$query);
  foreach($links as $lk){
    echo "|------ $lk ----|";
    $xml = simplexml_load_string($lk) or die("feed not loading");
    var_dump($xml);
  }
}
?>

```

There doesn't seem to be a `get()` function in PHP, you want `file_get_contents`. Also, your loops were nested incorrectly and you needed `var_dump($xml)` to dump the variable `$xml` not `var_dump(xml);` which tries to dump the string `xml`. You're not getting a file here but a string, so you need `simplexml_load_string`. Finally, you can't have unescaped whitespace in URLs.

Next time, the way to debug this sort of thing is to have it print the URL it is generating and try to visit it manually. You would have seen the spaces and noticed that the working example you show uses `%20` instead of space.
> 0 votes

# Answer
I don't know anything about php, but when I copy your code into here it isn't happy about the variable name `$query.`

This code appears to produce valid urls, but I can't make the `get` part work:

```
$eutils = "http://eutils.ncbi.nlm.nih.gov/entrez/eutils"; $db = "gene";$gene_symbols = array("A20","HBB","ABO");
$species = array("homo sapiens", "Escherichia coli","Saccharomyces cerevisiae"); $links = array();

foreach ($gene_symbols as $key => $GS){
   $SP = $species[$key];
   $SearchTerm = "$GS AND $SP [orgn] AND alive";

   $esearch = "$eutils/esearch.fcgi?db=$db&retmode=xml&retmax=10&sort=relevance&term=";

   echo $esearch.$SearchTerm;
}

http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=gene&retmode=xml&retmax=10&sort=relevance&term=A20 AND homo sapiens [orgn] AND alivehttp://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=gene&retmode=xml&retmax=10&sort=relevance&term=HBB AND Escherichia coli [orgn] AND alivehttp://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=gene&retmode=xml&retmax=10&sort=relevance&term=ABO AND Saccharomyces cerevisiae [orgn] AND alive

```
> 0 votes

# Answer
```
$eutils = "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi";

   $databank = "gene";
   $gene_symbols = array("A20","HBB","ABO");

   $species = array("homo sapiens", "Escherichia coli","Saccharomyces cerevisiae");

foreach($gene_symbols as $value => $GS){
    $Species = $species[$value];

    $esearch = "$GS%AND%alive[prop]%AND+$Species[orgn]";
    $esearch_test = "$eutils?db=$databank&retmode=xml&retmax=10&sort=relevance&term=$esearch";
    $result = file_get_contents($esearch_test);
    echo "Query: " . $esearch_test. "<br \>\n";
    echo "Result: " . $result . "<br \>\n";

    $xml = simplexml_load_string($result);
    echo "Count:<br \>\n";
    echo $xml->Count;
    $C =  $xml->Count;
    echo $C;
}

```
> 0 votes

---
Tags: phylogenetics
---


# Question
Title: What is cellranger doing in comparison to other methods?
I've recently started working with the 10X-Genomics platform with Illumina (MiSeq and HiSeq) for single-cell RNA-Seq. I've been recommended the "cellranger" (version 2.1.0) which I understand handles the barcoding of the platform and performs processing of the data. For example:

1) convert raw Illumina .bcl files into .fastq (demultiplex)

```
cellranger mkfastq

```

2) perform alignment, barcode counting, and generate gene-barcode-martices

```
cellranger count

```

However, it is unclear which underlying bioinformatics procedure (with which parameters) are being used here. For this pipeline, I wish to know how the trimming, alignment and counts are being performed. This is important to understand whether further normalisation and quality checking is necessary. What steps would I need to do to perform a similar analysis? If for example, I wish to document in more detail what has been performed here in a reproducible manner or reproduce the analysis pipeline on a different platform (such as Nadia) to compare them.

Ideally, I would like to know the analogous pipeline using open-source tools but at the moment I am first seeking to understand what is going on under the hood when I am running the `cellranger`.

# Answer
cellranger mkfastq is not necessary anymore. It used to be that the cellranger software wanted the reads to be interleaved, and you could use cellranger to do that for you if you couldn't do it yourself. Newer versions of cellranger will take the fastq files just like Illumina's bcl2fastq makes them.

Cellranger count aligns the reads, filters away duplicates based on UMIs, tries to figure out which cell barcodes really captured cells, and tells you for each cell barcode how many reads hit each gene. You can go through the .bams, and see which reads were counted for what gene and what barcode, based on the flags and other tags in the .bam

If you want another pipeline to compare to, you can use the McCarroll lab's DropSeq pipeline, which is a bunch of little scripts, around the output of STAR aligner (just like 10XGenomics uses). I have put the same data sets through both 10X and Dropseq, and they return about the same thing.
> 4 votes

# Answer
The best approach would be to contact 10x Genomics support: support@10xgenomics.com . They usually respond within a few hours.

Otherwise, you can check the methods from the first 10x Genomics paper:

> The Cell Ranger Single-Cell Software Suite was used to perform sample demultiplexing, barcode processing and single-cell 3′ gene counting (http://software.10xgenomics.com/single-cell/overview/welcome). First, sample demultiplexing was performed based on the 8 bp sample index read to generate FASTQs for the Read1 and Read2 paired-end reads, as well as the 14 bp GemCode barcode. Ten basepair UMI tags were extracted from Read2 (14 libraries were made with 5 bp UMI tags, as noted in Supplementary Table 1, due to an earlier iteration of the methods. For these samples, 5 bp UMI tags were extracted from Read2.). Then, Read1, which contains the cDNA insert, was aligned to an appropriate reference genome using STAR35. For mouse cells, mm10 was used and for human cells, hg19 was used. For samples with mouse and human cell mixtures, the union of hg19 and mm10 were used. For ERCC samples, ERCC reference (https://tools.thermofisher.com/content/sfs/manuals/cms\_095047.txt) was used.
>
> Next, GemCode barcodes and UMIs were filtered. All of the known listed of barcodes that are 1-Hamming-distance away from an observed barcode are considered. Then, the posterior probability that the observed barcode was produced by a sequencing error is computed, given the base qualities of the observed barcode and the prior probability of observing the candidate barcode (taken from the overall barcode count distribution). If the posterior probability for any candidate barcode is at least 0.975, then the barcode is corrected to the candidate barcode with the highest posterior probability. If all candidate sequences are equally probable, then the one appearing first by lexical order is picked.
>
> UMIs with sequencing quality score \>10 were considered valid if they were not homopolymers. Qual=10 implies 90% base call accuracy. A UMI that is 1-Hamming-distance away from another UMI (with more reads) for the same cell barcode and gene is corrected to the UMI with more reads. This approach is nearly identical to that in Jaitin et al.4, and is similar to that in Klein et al.8 (although Klein et al.8 also used UMIs to resolve multimapped reads, which was not implemented here).
>
> Last, PCR duplicates were marked if two sets of read pairs shared a barcode sequence, a UMI tag, and a gene ID (Ensembl GTFs GRCh37.82, ftp://ftp.ensembl.org/pub/grch37/release-84/gtf/homo\_sapiens/Homo\_sapiens.GRCh37.82.gtf.gz and GRCm38.84, ftp://ftp.ensembl.org/pub/release-84/gtf/mus\_musculus/Mus\_musculus.GRCm38.84.gtf.gz, were used). Only confidently mapped (MAPQ=255), non-PCR duplicates with valid barcodes and UMIs were used to generate gene-barcode matrix.
>
> Cell barcodes were determined based on distribution of UMI counts. All top barcodes within the same order of magnitude (\>10% of the top nth barcode, where n is 1% of the expected recovered cell count) were considered cell barcodes. Number of reads that provide meaningful information is calculated as the product of four metrics: (1) valid barcodes; (2) valid UMI; (3) associated with a cell barcode; and (4) confidently mapped to exons.

Of course, they have since updated the library prep and the Cell Ranger software, so the exact protocol is somewhat different.
> 3 votes

---
Tags: rna-seq, sequence-analysis, single-cell, 10x-genomics, cellranger
---


# Question
Title: How to add OVERLAP_BP field to VEP output?
I'm using the Ensembl Variant Effect Predictor (VEP) to annotate VCF files, and I am trying to add some additional fields to the output, among these is `OVERLAP_BP`, which is listed here in available fields:

> OVERLAP\_BP - Number of base pairs overlapping with the corresponding variation feature

It's clear that some (many? all?) of these additional fields are only enabled by adding extra options to the VEP command (as listed on the RHS of the options table).

For example, the `ZYG` field is made available by adding the `--individual` flag to the VEP call.

But `OVERLAP_BP` doesn’t appear in the output fields column on that page, and when I add it to the output using `--fields`, it is consistently empty.

What options are required to make the `OVERLAP_BP` output field available to me in VEP?

# Answer
> 1 votes
I also opened this as a Github issue, and received the following response:

> It seems there are some inconsistencies on our end. I'm looking into it. Until the fix, using `--fields OverlapBP,OverlapPC` works with any output format, while `--fields OVERLAP_BP,OVERLAP_PC` works only with default output format.
>
> A small note that this overlap data is available for structural variants only.

---
Tags: vcf
---


# Question
Title: Manually define clusters in Seurat and determine marker genes
I want to define two clusters of cells in my dataset and find marker genes that are specific to one and the other. Is there a way to do this in `Seurat`? Say, if I produce two subsets by the `SubsetData` function, is there a way to feed them into some other function that would calculate marker genes? If not, what other packages would you recommend for doing that?

If you look here:

https://satijalab.org/seurat/seurat\_clustering\_tutorial\_part2.html

I just need a way to define `ident` myself, the number of levels (2) and assign numbers to each cell (0, 1), and then run `DE` between `0` and `1` clusters which is obvious how to do afterwards.

# Answer
> 3 votes
I think you are looking to `FindAllMarkers` function from `Seurat`. As you said, you just have to define your `ident`, that have to have the structure of a table (cell names as names and cluster as value):

```
pident=as.factor(clusters)
names(pident)=cellNames
object1@ident=pident

```

And then run the `FindAllMarkers` function:

```
FindAllMarkers(object1, min.pct = 0.25, min.diff.pct = 0.25)

```

You can specify several parameters in this function (type of DE to perform, thresholds of expression, etc).

# Answer
> 3 votes
Seurat has functions for adding metadata and setting identities. Get unique cell names:

```
cell.labels <- seuratobject@ident

```

Replace column and its name with your cluster labels (e.g.), then:

```
seuratobject <- AddMetaData(seuratobject, metadata=cell.labels)
seuratobject <- SetAllIdent(seuratobject, id='yourclusterlabels')

```

Because you want to contrast two clusters against each other, I suggest using `FindMarkers()` as opposed to `FindAllMarkers()`:

```
FindMarkers(object, ident.1, ident.2)

```

It can also compare combinations of clusters.

---
Tags: r, single-cell, seurat
---


# Question
Title: Is there a standard definition for "assembly polishing"?
1. Is there a standard definition for "assembly polishing" in the field?
2. Is there a standard definition for what polishing algorithms do?

My understanding of "polishing" is strongly influenced by Pilon:

Whereby users input a reference FASTA along with the BAM of aligned reads against this reference FASTA. The polishing algorithm then does alignment to discover local misassemblies and other alignment inconsistencies. These are then "corrected", i.e. "polished".

# Answer
> 6 votes
In long-read assembly, "polish" refers to the step to improve the base accuracy of contig sequences. I believe the terminology was originated from the HGAP paper:

> The final consensus–calling algorithm Quiver, which takes into account all of the underlying data and the raw quality values inherent to SMRT sequencing, then ***polishes*** the assembly for final consensus accuracies in excess of 99.999% (QV of \>50)...

Pilon was published more than a year later and in fact, neither its paper nor README mentions word "polish".

By the way, "error correction" refers to the step to fix sequencing errors in raw reads, not in contigs. In Canu/HGAP/Falcon, that step is done by pbdagcon or falcon\_sense.

# Answer
> 5 votes
**Disclaimer**, this answer is based on feeling I got from talks and papers, but I do not have any hard reference supporting it.

I believe that genome polishing is a technique that was introduced for correction of individual bases of long noisy read assemblies - i.e. fixing SNPs and short indels. Guessing a correct base from long reads during assembly seems to be bit more complicated task in comparison to short read assemblies, therefore it is isolated into separated algorithms like Quiver/Arrow, Nanopolish (for correction by long reads) or the tool you mentioned - Pilon (correction by short reads).

Miss-assemblies could be corrected even prior long reads for instance by matching the assembly to optical mapping, or by tools like REAPR or TGnet. However, I have never heard term "genome polishing" used for this type of genome correction.

---
Tags: sequence-alignment, assembly
---


title: Protein sequence from patient data

text:
# Question
Title: Protein sequence from patient data
Currently, I am working on NGS data and my aim is to get significance prediction of variants present in the vcf file. As we know about SIFT Score for significance score prediction, I am trying to understand how this score works. When I read its reference paper, in the method section, the procedure to determine the sift score starts with patient protein query sequence.

My question is how to get patient protein sequence from patient data (either vcf file or bam file) (which is not explained in the paper). Rest I can follow by other algorithms like PSI-BLAST & MOTIF.

Any suggestions....

Thanks

# Answer
> 5 votes
Ensembl now have a tool called Haplosaurus (still in beta) which will convert your phased VCF into actual protein sequences.

> ## Haplosaurus
>
> `haplo` is a local tool implementation of the same functionality that powers the Ensembl transcript haplotypes view. It takes phased genotypes from a VCF and constructs a pair of haplotype sequences for each overlapped transcript; these sequences are also translated into predicted protein haplotype sequences. Each variant haplotype sequence is aligned and compared to the reference, and an HGVS-like name is constructed representing its differences to the reference.
>
> This approach offers an advantage over VEP's analysis, which treats each input variant independently. By considering the combined change contributed by all the variant alleles across a transcript, the compound effects the variants may have are correctly accounted for.
>
> `haplo` shares much of the same command line functionality with `vep`, and can use VEP caches, Ensembl databases, GFF and GTF files as sources of transcript data; all `vep` command line flags relating to this functionality work the same with `haplo`.
>
> ## Usage
>
> Input data must be a VCF containing phased genotype data for at least one individual; no other formats are currently supported.
>
> When using a VEP cache as the source of transcript annotation, the first time you run haplo with a particular cache it will spend some time scanning transcript locations in the cache.
>
> ```
> ./haplo -i input.vcf -o out.txt -cache
>
> ```

# Answer
> 3 votes
You can use SIFT in two different ways. One way is manually via their website. This is done by downloading the protein sequence (fasta) from a reliable repository (e.g., uniprot or NCBI), and change the amino acid of interest. Then calculate the SIFT score.

The other way, which is probably what you are interested in is with NGS data. Here normally SIFT is already known (precalculated) for all coding SNPs, and you can use that info directly during annotation of your variants. I don't know which variant annotation pipeline you are using, but for example annovar has SIFT data available at their databases repository.

---
Tags: vcf, ngs, variation
---


title: Large dataset normalization for PCA

text:
# Question
Title: Large dataset normalization for PCA
I need to normalize a large (400Mb) dataset for doing PCA analysis. I want to use `scran` for doing that:

```
biocLite("SingleCellExperiment")
biocLite("scran")
library(SingleCellExperiment)
library(scran)
list_of_sce <- list()
# Looping though the UMI_count 'split_factor' columns at a time
split_factor = 500
for(i in seq(1,ncols, split_factor)) {
  num_loop = floor(i / split_factor) + 1
  idx = ncols
  if (i + split_factor < ncols) {
  idx = i + split_factor
  }
  sce <- SingleCellExperiment(list(counts=UMI_count[, i : idx]))
  # Normalization of the dataset containing
  # heterogenous cell data (different cell types)
  clusters <- quickCluster(sce) ## <- ** error appears here **
  sce <- computeSumFactors(sce, cluster=clusters)
  list_of_sce[[num_loop]] <- sce
}
cbind(list_of_sce)

```

However, when I try to run this it produces an error at the `quickCluster` step:

```
Error: cannot allocate vector of size 23.7 Gb

```

It feels like `cbind` at the end is not the right way to go, and I need to somehow combine `sce` and re-normalize them again somehow. The question is how to do that in a correct way: normalizing large dataset for PCA analysis when it can not be fully loaded into memory?

# Answer
I normalized a high throughput dataset for a school project using `DESeq` library using the script bellow. The code is based on a lesson I had. My goal was determine the over expressed genes, but the normalization step should be the same.

```
#####################################################################
### DATA PRE-PROCESSING: Normalising gene expression distributions
#####################################################################
# Genes that are not expressed at a biologically meaningful level
# in any condition should be discarded to reduce the subset of genes
# to those that are of interest, and to reduce the number of tests carried out
# downstream when looking at differential expression.
# Using a nominal CPM value of 1 (which is equivalent to a log-CPM value of 0
# a gene is deemed to be expressed in a given sample if its transformed count
# is above this threshold, and unexpressed otherwise.
# Genes must be expressed in at least one group (or roughly any three samples
# across the entire experiment) to be kept for downstream analysis.

library("DESeq")

keep.exprs <- rowSums(cpm>1)>=3 # is filtering is ok
DGEfiltered <- DGEobject[keep.exprs,, keep.lib.sizes = FALSE]

## Normalisation generally required to ensure that expression distributions
## of each sample are similar across the entire experiment.
## Any plot showing the per sample expression distributions, such as a density or boxplot,
## is useful in determining whether any samples are dissimilar to others.
## Normalisation by the method of trimmed mean of M-values (TMM)
## is performed using the calcNormFactors function in edgeR.
# cpmF = filtered subset # log.cpmF = log cpm filtered subset
# Copy unnormalized, filtered data to show in boxplot
DGEf.notnorm <- DGEfiltered # normaliseerde
DGEf.norm <- calcNormFactors(DGEfiltered, method = "TMM")
DGEf.norm$samples$norm.factors

# --> scaling factors relatively close to 1 --> mild TMM normalisation

## Boxplots of log-CPM values showing expression distributions
## for unnormalised data (A) and normalised data (B) for each sample
par(mfrow = c(1,1))
par(mar = c(7,4,2,1))
# Boxplot unnormalised data
DGEf.notnorm$samples$norm.factors
# [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1
log.cpmF.notnorm <- cpm(DGEf.notnorm, log = TRUE)
boxplot(log.cpmF.notnorm, las = 2, col = sample.col, cex = 0.9,
        main = "A. Unnormalised data", ylab = "Log-cpm")
# Boxplot normalised data
DGEf.norm$samples$norm.factors

log.cpmF.norm <- cpm(DGEf.norm, log = TRUE)
boxplot(log.cpmF.norm, las = 2, col = sample.col, cex = 0.9,
        main = "B. Normalised data", ylab = "Log-cpm")
# --> only when scaling factors are not close to 1 you will see a big difference

```
> 3 votes

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


# Question
Title: NCBI xml-format Printing nodeValues from tags with the same names
I would like to retrieve all Other-source\_anchor on the level of DNA binding from within a Gene Ontology XML file using php.

If I try this with my code, I only get Gene Ontology back and not the nodeValues in Other-source\_anchor. Could someone tell me how to do this?

# Xml file were I wwant to get the nodeValues from

```
<Gene-commentary_heading>GeneOntology</Gene-commentary_heading>
      <Gene-commentary_source>
        <Other-source>
          <Other-source_pre-text>Provided by</Other-source_pre-text>
          <Other-source_anchor>GOA</Other-source_anchor>
          <Other-source_url>http://www.ebi.ac.uk/GOA/</Other-source_url>
        </Other-source>
      </Gene-commentary_source>
      <Gene-commentary_comment>
        <Gene-commentary>
          <Gene-commentary_type value="comment">254</Gene-commentary_type>
          <Gene-commentary_label>Function</Gene-commentary_label>
          <Gene-commentary_comment>
            <Gene-commentary>
              <Gene-commentary_type value="comment">254</Gene-commentary_type>
              <Gene-commentary_source>
                <Other-source>
                  <Other-source_src>
                    <Dbtag>
                      <Dbtag_db>GO</Dbtag_db>
                      <Dbtag_tag>
                        <Object-id>
                          <Object-id_id>3677</Object-id_id>
                        </Object-id>
                      </Dbtag_tag>
                    </Dbtag>
                  </Other-source_src>
                  <Other-source_anchor>DNA binding</Other-source_anchor>

```

# code:

```
$esearch_test = "https://www.ncbi.nlm.nih.gov/gene/7128?report=xml&format=text";
$result = file_get_contents($esearch_test);
$xml = simplexml_load_string($result);
$doc = new DOMDocument();
$doc = DOMDocument::loadXML($xml);
$c = 1;
foreach($doc->getElementsByTagName('Gene-commentary_heading') as $node) {
    if ($node->textContent =="GeneOntology"){
        echo "<pre>"."$c: ".$node->nodeName."</pre>";
        echo "<pre>"."$c: ".$node->nodeValue."</pre>";
        $cn = 1;
        foreach ($node->getElementsByTagName('Other-source_anchor') as $subnode){
            echo "<pre>"."$c: ".$subnode->nodeName."</pre>";
            echo "<pre>"."$c: ".$subnode->nodeValue."</pre>";

            echo $cn++;
        }
    }
    $c++;
}

```

# Answer
use a XSLT processor (http://php.net/manual/en/class.xsltprocessor.php) with the following stylesheet.

```
<?xml version='1.0' encoding="UTF-8"?>
<xsl:stylesheet xmlns:xsl='http://www.w3.org/1999/XSL/Transform' version='1.0'>
    <xsl:output method="html"/>

    <xsl:template match="/">
        <div>
            <xsl:apply-templates select="//Gene-commentary_heading[ text() = 'GeneOntology']"/>
        </div>
    </xsl:template>

    <xsl:template match="Gene-commentary_heading">
        <pre><xsl:value-of select="name(.)"/></pre>
        <pre><xsl:value-of select="."/></pre>
        <xsl:apply-templates select="..//Other-source_anchor"/>
    </xsl:template>

    <xsl:template match="Other-source_anchor">
        <pre><xsl:value-of select="name(.)"/></pre>
        <pre><xsl:value-of select="."/></pre>
    </xsl:template>

</xsl:stylesheet>

```
> 3 votes

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