# Question
Title: What process and input data is required for a cellranger reference transcriptome?
I'm analysing single-cell RNA-Seq data using the 10X Genomics `cellranger` platform. While they provide reference data for Mouse and Humans, other species require a reference "transcriptome" to be generated with the cellranger `mkref` tool. However, I've found it to be very temperamental giving errors such as:

```
gffread genome_GFF3_genes.gff -T -o genome_genes.gtf

cellranger mkgtf genome_genes.gtf genome_genes.filtered.gtf \
                   --attribute=gene_biotype:protein_coding \
                   --attribute=gene_biotype:lincRNA \
                   --attribute=gene_biotype:antisense \
                   --attribute=gene_biotype:IG_LV_gene \
                   --attribute=gene_biotype:IG_V_gene \
                   --attribute=gene_biotype:IG_V_pseudogene \
                   --attribute=gene_biotype:IG_D_gene \
                   --attribute=gene_biotype:IG_J_gene \
                   --attribute=gene_biotype:IG_J_pseudogene \
                   --attribute=gene_biotype:IG_C_gene \
                   --attribute=gene_biotype:IG_C_pseudogene \
                   --attribute=gene_biotype:TR_V_gene \
                   --attribute=gene_biotype:TR_V_pseudogene \
                   --attribute=gene_biotype:TR_D_gene \
                   --attribute=gene_biotype:TR_J_gene \
                   --attribute=gene_biotype:TR_J_pseudogene \
                   --attribute=gene_biotype:TR_C_gene

cellranger mkref --genome="reference-name" \
                 --fasta="transcriptome.fasta" \
                 --genes="genome_genes.filtered.gtf"

```

> Fatal INPUT FILE error, no valid exon lines in the GTF file: reference\_name/genes/genes.gtf Solution: check the formatting of the GTF file. Most likely cause is the difference in chromosome naming between GTF and FASTA file.
>
> subprocess.CalledProcessError: Command '\['STAR', '--runMode', 'genomeGenerate', '--genomeDir', 'reference\_name/star', '--runThreadN', '1', '--genomeFastaFiles', 'reference\_name/fasta/genome.fa', '--sjdbGTFfile', 'reference\_name/genes/genes.gtf', '--limitGenomeGenerateRAM', '17179869184', '--genomeSAsparseD', '1', '--genomeSAindexNbases', '12', '--genomeChrBinNbits', '10'\]' returned non-zero exit status 104
>
> subprocess.CalledProcessError: Command '\['STAR', '--runMode', 'genomeGenerate', '--genomeDir', 'reference\_name/star', '--runThreadN', '1', '--genomeFastaFiles', 'reference\_name/fasta/genome.fa', '--sjdbGTFfile', 'reference\_name/genes/genes.gtf', '--limitGenomeGenerateRAM', '17179869184', '--genomeSAsparseD', '1', '--genomeSAindexNbases', '12', '--genomeChrBinNbits', '10'\]' returned non-zero exit status -11

What input data is required for the `cellranger` to construct a reference? What format do the fasta and gtf files need be in? Is it possible to use a "transcriptome" rather than a "genome" as a reference here? For example, I have a transcriptome reference with fasta names `>name` for each transcript rather than each chromosome.

# Answer
Your problem is caused by using the transcriptome fasta file rather than the genome fasta file. You've already given it transcriptome information with `genome_genes.filtered.gtf`, it needs the genomic sequence to go along with that.

As an aside, your GTF filtering command is unlikely to be useful unless you have a Gencode or Ensembl GTF for human. Have a look at the `gene_biotype`s in your file and decide for yourself what (if anything) needs to be filtered. For example, I made a cellranger reference last night for dm3 (the old fruit fly genome) with:

```
cellranger mkgtf annotation.gtf genes.gtf \
    --attribute=gene_biotype:protein_coding \
    --attribute=gene_biotype:ncRNA \
    --attribute=gene_biotype:snRNA \
    --attribute=gene_biotype:snoRNA \
    --attribute=gene_biotype:pre_miRNA

```

That should be customized for relevance. It's quite possible that your GTF doesn't even have `gene_biotype` attributes, that's specific to Ensembl I think.
> 4 votes

---
Tags: scrnaseq, reference-genome, transcriptome, star, cellranger
---


# Question
Title: Automating download of results of Drug and Drug Target Mapping
The drug and drug target mapping tool shows the **top 3** PDB IDs for each drug based on the drug target sequence similarity search and ligand id search. Is there a way to get all the PDB ids of the drug-target PDB files and not only the top 3?

I know you can get this information after clicking the "find" button under the last column of the output table (the column named "All PDB"), but I have a query of several hundreds of drugs and I'm not willing to do it manually for each result.

# Answer
Simpler option would be to not use the query interface and directly parse the CSV file found on that page.

I just read @llrs comment and I guess it's very much the same spirit as my own - download the data and query it programmatically. I see now (was on mobile earlier) that the 'FIND' URL links out to a BLAST search. so it might be as easy as retrieving the relevant protein sequence and replicating the URL. Or if you want to go a different route and finding all PDB identifiers with the given ligand, cross-reference the above data with what you'd find here
> 0 votes

---
Tags: proteins, pdb
---


# Question
Title: how to do analysis of a table content SNP?
I have a table with information about SNPs, the table has a total of 55127 rows.

How can I know the number of gene mutations in the table and if a gene is present?

The table is in xls format.

# Answer
For working with a relatively small dataset in tabular format, I would recommend dplyr. You can read in a `excel` file using readxl

```
library(readxl)
library(dplyr)

df <- read_excel(FILENAME)

```

To get the number of mutations per gene:

```
df %>% group_by(`Gene name`) %>%
  summarise(num_mutations=n())

```

If you want to know if there are mutations in a particular gene:

```
df %>% filter(`Gene name` == QUERY)

```
> 1 votes

---
Tags: r, snp, mutations
---


# Question
Title: Transform a R data.frame into a krona plot without krona tools?
KronaTools lacks a R pipeline to convert any data.frame into a krona plot. There is a way to go via a phyloseq-class object, but coming from a data.frame I did not find any documentation on how to transform these.

# Answer
First download any Krona plot (.html) file or make one by using the Excel sheet from the Krona Package. Open in text-editor and remove anything below "datasets" tag and/or the first "node" tag till the end of the document.

Here is the path to go: 1) get your data.frame in hierachical order 2) transform to data.tree object 3) transform to xml file 4) gsub until all the tags match krona annotation

```
options(scipen=999)
options(stringsAsFactors=FALSE)
library(data.tree)
library(xml2)
library(magrittr)

#get your data into a hieracial format in a data.frame
t <- structure(list(freq = c(34L, 127L, 42L, 123L, 343L, 13L, 104L,
21L, 135L, 91L, 167L, 148L, 90L, 46L, 25L, 79L, 1726L), king.n = c("",
"Fungi", "Fungi", "Fungi", "Fungi", "Fungi", "Fungi", "Plantae",
"Plantae", "Plantae", "Plantae", "Plantae", "Plantae", "Plantae",
"Plantae", "Plantae", "Plantae"), class.n = c("", "Dothideomycetes",
"Dothideomycetes", "Dothideomycetes", "Dothideomycetes",     "Dothideomycetes",
"Dothideomycetes", "Liliopsida", "Magnoliopsida", "Magnoliopsida",
"Magnoliopsida", "Magnoliopsida", "Magnoliopsida", "Magnoliopsida",
"Magnoliopsida", "Magnoliopsida", "Magnoliopsida"), fam.n = c("",
"", "Cladosporiaceae", "Didymellaceae", "Pleosporaceae", "Pleosporaceae",
"Pleosporaceae", "Poaceae", "", "Asteraceae", "Asteraceae", "Asteraceae",
"Asteraceae", "Fabaceae", "Hydrophyllaceae", "Hypericaceae",
"Malvaceae"), gen.n = c("", "", "Cladosporium", "", "Alternaria",
"Pyrenophora", "Pyrenophora", "Phleum", "", "", "Achillea", "Centaurea",
"Jurinea", "Trifolium", "Phacelia", "Hypericum", "Malva"), spec.n = c("",
"", "cladosporioides", "", "", "", "tritici-repentis", "pratense",
"", "", "", "cyanus", "", "", "", "", "")), .Names = c("freq",
"king.n", "class.n", "fam.n", "gen.n", "spec.n"), row.names = c(1L,
3L, 4L, 5L, 8L, 9L, 11L, 15L, 16L, 17L, 18L, 19L, 21L, 22L, 25L,
26L, 28L), class = "data.frame")

#if any sequences can not be assigned to the highest hieracial level give them a random name like "NA"
t[t[,"king.n"]=="","king.n"] <- "NA"
#define the hieracy (do not include the counting var here), the first space is the root category
t$pathString <- paste(" ",t$king.n,t$class.n,t$fam.n,t$gen.n,t$spec.n,sep = "/")
#transform to a data.tree
t <- as.Node(t)
#define the sequence frequency calculations
t$Do(function(x) x$freq <- Aggregate(x, "freq", sum), traversal = "post-order")
print(t,"freq")
#get all frequencies for each node
repval <- t$Get("freq")

#tranform into a nested list and afterwards into a xml document to get the nesting
t <- t %>% as.list(.,mode="simple",unname=T)
t <- as_xml_document(list(root=t))%>%as.character

#from here it is just a bunch of regex
#split the xml file into one line for each <
temp <- strsplit(t,"\n")[[1]][2] %>% strsplit(.,"(?=<)",perl=T) %>% .[[1]]
temp <- paste0(temp[c(T,F)],temp[c(F,T)]) %>% gsub("[0-9]","",.)
#remove the names nodes and replace by krona annotation
#replace the closing named tags by generic closing "node" tags
temp <- gsub("</(.*)?>","</node>",temp) %>% gsub(">.*$",">",.)
#replace the named opening tags by generic "node" tags with the "name" attribute, replace all (because they are incomplete) numeric values with a placeholder
temp <- gsub("<(?=[^/])(.*)?>",
     '<node name="\\1"><magnitude><val>__\\1</val></magnitude><score><val>1</val></score>',temp,perl=T)
temp %<>% gsub("root"," ",.)
#find the placeholders to insert the frequencies and execute
repos <- grep("__",temp)
for (i in repos)
  temp[i] <- gsub("^(.*)(__.*?)(?=<)(.*)",paste0("\\1",repval[which(repos==i)],"\\3"),temp[i],perl=T)

#add the end of the .html file
temp <- c(temp,"</krona></div>","</body>","</html>")  %>% unlist
tout <- paste(temp,collapse="\n")

#append your data to the Krona plot
write(tout,file="Krona1.html",append=TRUE)

```

Et voilá!
> 3 votes

---
Tags: r
---


title: GFF file too big to load into R

text:
# Question
Title: GFF file too big to load into R
I want to load a GFF (annotation) file into an R data frame in order to extract some information from it, e.g. the location of transposons in the genome.

I used the following code:

```
library(rtracklayer)
gff <- rtracklayer::import('c_elegans.PRJNA13758.WS263.annotations.gff3.gz')
Error: cannot allocate vector of size 210.2 Mb

```

It seems that the GFF3 file is too big for loading into R. Do you know if this is limited by R or by my computer's memory?

I ran this on a GNU/Linux computer.

# Answer
It depends of course on your available RAM on your computer.

Under Windows you can increase the available RAM in R with memory limit:

```
> memory.limit()
[1] 8070
> memory.limit(size=10000)
[1] 10000
> memory.limit()
[1] 10000

```

But that is not necessary on linux/unix.
> 3 votes

---
Tags: r, bioconductor, sequence-annotation, errors
---


# Question
Title: single stranded cDNA vs double-stranded cDNA
Single stranded cDNA is used by race pcr and rt-pcr while next-gen sequencing (Illumina) requires double stranded cDNA. What is the reason behind this ?

# Answer
Sequencing also requires single stranded DNA, otherwise the sequencing primers can’t bind. Sure there are steps in the protocol where double stranded material exists, but that’s transient.
> 1 votes

---
Tags: phylogenetics
---


# Question
Title: How to compare the contents of a column of the same data table
I have this table:

and I want to get the rows that are equal to the first three columns, like this:

I've tried these functions, but when I get the index of the lines, r doesn't give the output that I want:

```
df$obj<-sapply(c("sample1", "sample2", "sample3"), function (comparison) {
    df$sample1 == df[comparison]
})

```

# Answer
Loop through 3 columns by row using *apply*, then convert *toupper* case (as we want `a` same as `A`), get *unique*, and get *lengths*, if length is 1 then all 3 values are equal:

```
df1[ lengths(
  apply(df1[, 1:3], 1, function(i) unique(toupper(i)))
  ) == 1, ]

```

This would scale better if we want to compare 4, 10, n(?) columns.
> 4 votes

# Answer
No need to use *dplyr*, we can do this using *base*, see below example:

```
# reproducible example data
df1 <- data.frame(
  t(data.frame(
    one = c("A", "A", "A", "C", "C", "C"),
    two = c("A", "A", "A", "A", "A", "C"),
    thr = c("A", "A", "A", "A", "A", "A"),
    fou = c("C", "A", "A", "A", "A", "A"))),
  stringsAsFactors = FALSE)

df1
#     X1 X2 X3 X4 X5 X6
# one  A  A  A  C  C  C
# two  A  A  A  A  A  C
# thr  A  A  A  A  A  A
# fou  C  A  A  A  A  A

# solution, subset using logical index:
selection <- df1[df1[, 1] == df1[, 2] & df1[, 1] == df1[, 3], ]

selection
#     X1 X2 X3 X4 X5 X6
# one  A  A  A  C  C  C
# two  A  A  A  A  A  C
# thr  A  A  A  A  A  A

```
> 3 votes

---
Tags: r
---


# Question
Title: Why isn't there a standard method to convert GFA to JSON?
I'm confused why the GFA format doesn't have a standard JSON parser. It appears to be related to a decision made several years ago, here.

The objection appears to be "JSON isn't necessary, as it's just a simple linked list. Users could write their own". However, I would think this is a useful feature.

There are applications I could think of whereby a JSON would be necessary, or at least very convenient.

(1) Is there a standard way to convert GFA into JSON? Or do users write their own scripts?

(2) Is there a conceptual reason why GFA doesn't convert into JSON?

# Answer
Short answer: no, there isn't a "standard" method, and there won't be.

Long answer: like TAB-delimited and XML, JSON alone is not a specific format. You have to define a schema to give meanings to data. However, unlike TAB-delimited "format" and XML, there is not an official way to define JSON schema. There have been various attempts but all of them are complicated and none of them are widely adopted. As you can't write down a JSON schema for GFA in the first place, there won't be a "standard" way to convert GFA to JSON. The best you can do is to roll out your own in-house solution, unfortunately.
> 4 votes

---
Tags: assembly
---


# Question
Title: Biopython Phylogenetic Tree replace branch tip labels by sequence logos
Having recently constructed a lot of phylogenetic trees with the module TreeConstruction from Phylo package from Biopython, I've been asked to replace the branch tip labels by the corresponding sequence logos (which I have in the same folder). I thought that it would be more efficient to make a code to generate the logo-trees automatically, as I would have to make a lot of them.

The first idea that I came up with was to see whether the functions used to build the tree had an argument to replace the branch tip labels or to remove them, which I couldn't find. Therefore, I removed the branch tip labels by setting their font size to 0: (following is the code to build the tree)

```
# Modules to build the tree
from Bio.Phylo.TreeConstruction import DistanceCalculator, DistanceTreeConstructor
from Bio.Phylo import draw
from Bio import Phylo, AlignIO
import subprocess
import matplotlib
import matplotlib.pyplot as plt

alignment = AlignIO.read('MotifSeqAligned.fasta', 'fasta') # reading the alignment file

calculator = DistanceCalculator('ident')
dm = calculator.get_distance(alignment) # distance matrix

constructor = DistanceTreeConstructor()
tree = constructor.nj(dm) # build with neighbour joining algorithm a tree from dm

Phylo.write(tree, 'TreeToCutOff.nwk', 'newick')

plt.rc('font', size=0)          # controls default text sizes #HERE IS THE SETTING FOR THAT ALLOWS ME TO HIDE THE BRANCH TIP LABELS
plt.rc('axes', titlesize=14)     # fontsize of the axes title
plt.rc('xtick', labelsize=10)    # fontsize of the tick labels
plt.rc('ytick', labelsize=10)    # fontsize of the tick labels
plt.rc('figure', titlesize=18)   # fontsize of the figure title

draw(tree, do_show=False)
plt.savefig("TreeToCutOff.svg", format='svg', dpi=1200)

```

From this code I could get the tree:

Since I don't know how to get the y coordinates of the branches to add the logos one by one, I built a column of logos with matplotlib, that I intended then to paste on the tree in python. The code to build the column of logos is the following:

```
#Extract filename from newick
newickFile = open("TreeToCutOff.nwk", 'r').read()
orderedLogos = ["{}.eps".format(i) for i in re.split('(\W)', newickFile) if "Profile" in i]

#Initialize the figure
fig = plt.figure()

#Add each image one after the other in the right order
for i, files in enumerate(orderedLogos):
    img1 = mpimg.imread(files)
    ax1 = fig.add_subplot(len(orderedLogos), 1, 1+i)
    ax1.imshow(img1)
    ax1.set_xticks([])
    ax1.set_yticks([])

# plt.show()
plt.savefig("RowsOfLogos.svg", format='svg', dpi=1200)
plt.clf()
plt.cla()

```

Having my tree and the column of logos in .svg or .png, I couldn't find any way to stack them properly. My first idea was to use the library svgutils which seemed to be easy to handle, with the following code: (taken from svgutils tutorials)

```
import svgutils.transform as sg
# Assemble
#create new SVG figure
fig = sg.SVGFigure("14cm", "14cm")

# load figures
fig1 = sg.fromfile('TreeToCutOff.svg')
fig2 = sg.fromfile('RowsOfLogos.svg')

# get the plot objects
plot1 = fig1.getroot()
plot2 = fig2.getroot()
plot2.moveto(280, 100, scale=0.05)

# append plots and labels to figure
fig.append([plot1, plot2])

```

But the issue with the output was that the background of the column of logos was white and thus, I was pasting a huge white image with a thin column of logos on the tree. And I couldn't find a way to crop the column of logos with svgutils. I tried the module Image from PIL package to build a tree of logos from .png files but couldn't see the tree used as background after pasting the column of logos.

There may be a way to do what I'm aiming for with matplotlib (which would be to stack 2 files .png, and place the logos all the time at the same distance), but I couldn't work it out.

Does anyone know what the best solution is to make a tree of logos (as in the following image which I could only build manually with inkscape) with python libraries, allowing to automate the process without having to adapt the code depending on the number of branches?

Following is a subset of "MotifSeqAligned.fasta" containing the aligned sequences used to build the trees:

```
>ProfileCluster0.meme
---SSNDTTTCCAGGAAD-
>ProfileCluster1.meme
YBNRD---TTCYYGGAAT-
>ProfileCluster10.meme
---VDKDWTTCTYGGAAT-

```

With the 3 corresponding logos:

The full "MotifSeqAligned.fasta" and all the logos.eps (as I used them instead of .png which is the format asked by the forum) can be found here.

# Answer
To expand on my comment from yesterday. You could do this with the ETE Toolkit (I just copied one logo file rather than converting all 26 to `png`):

```
from ete3 import Tree, TreeStyle, faces
def mylayout(node):
    if node.is_leaf():
        logo_face = faces.ImgFace(str.split(node.name, '.')[0] + ".png") # this doesn't seem to work with eps files. You could try other formats
        faces.add_face_to_node(logo_face, node, column=0)
    node.img_style["size"] = 0  # remove blue dots from nodes

t = Tree("tree.nwk", format=3)
ts = TreeStyle()
ts.layout_fn = mylayout
ts.show_leaf_name=False # remove sequence labels
ts.scale =  10000 # rescale branch lengths so they are longer than the width of the logos
t.render("formatted.png", tree_style = ts, h=3000, w=3000) # you may need to fiddle with dimensions and scaling to get the look you want

```

If you want all of the logos lined up in a column add `aligned=True` to `faces.add_face_to_node`
> 8 votes

---
Tags: python, motifs, phylogenetics
---


# Question
Title: Compute copy number from cases and controls
I have some data on Copy Number Variation (SNP chip) for a population of samples.

In particular, I have a set of samples (considered as cases) which display a specific disease phenotype, and another set (considered as controls) which do not. The cases have not a matched-control. All the controls are taken as a random sample from the population that does not present the disease.

For both cases and controls I have the copy number of some regions. How can I compute the copy number change in cases, such as in GDC?

Should I use (as described in the previously linked page):

> performing tangent normalization, which subtracts variation that is found in a set of normal samples

and do you know of any tool to perform this computation?

The data I have is formatted in this way:

```
Sample_ID   chrom      start           end  CN
Sample11       19   11991477      12133823   1
Sample11        2   52260564      52431658   1
Sample12        7    5721757       5896192   3
Sample13       10    2269963       2473585   3

```

# Answer
GISTIC does exactly this, and is very possibly the tool used in your link above.

The input files need to be in a different format than what you currently have, but if you're using a SNP6 array, there are guides out there that will tell you how to get the proper files from your `.cel` files.
> 3 votes

---
Tags: cnv, copy-number, snp-chip
---


# Question
Title: How can I get a list of papers where a particular mutation has been mentioned?
I have recently tried a few tools (MutationFinder, tmVar, PubTator etc.) that extract mutations from a given text. However, for my work I need something that will take a list of mutations as an input and return a list of pubmed id in which the mutations has been found.

Example input: mut\_list = \[Y702F, L483Y\]

Example output: pub\_list = \[12215417, 18537569\]

Can anyone give any idea?

Thanks!

# Answer
> 5 votes
You could use the pubmed API to query it directly. Typically I use Biopython for my pubmed queries. There are some good examples in the cookbook and on other SO posts.

This example from the cookbook link above has the guts of what you need. You will just have to iterate over your list. Depending on the number of queries you want to do you might have issues getting throttled by NCBI

```
In this example, we will query PubMed for all articles having to do with orchids (see section 2.3 for our motivation). We first check how many of such articles there are:

>>> from Bio import Entrez
>>> Entrez.email = "A.N.Other@example.com"     # Always tell NCBI who you are
>>> handle = Entrez.egquery(term="orchid")
>>> record = Entrez.read(handle)
>>> for row in record["eGQueryResult"]: ...
        if row["DbName"]=="pubmed": ...
            print(row["Count"]) 463

Now we use the Bio.Entrez.efetch function to download the PubMed IDs of these 463 articles:

>>> from Bio import Entrez
>>> handle = Entrez.esearch(db="pubmed", term="orchid", retmax=463)
>>> record = Entrez.read(handle)
>>> handle.close()
>>> idlist = record["IdList"]

This returns a Python list containing all of the PubMed IDs of articles related to orchids

```

---
Tags: mutations
---


# Question
Title: Smallest group size for differential expression in limma (bulk RNA-Seq)
I am reading Smyth et al. (ref. 1). I want to run differential expression analysis on a bulk RNA-Seq dataset in which each group is composed by 2 samples. In the paper previously cited it is written that:

> Genes must be expressed in at least one group (or in at least three samples across the entire experiment, where **three was chosen as this is the smallest group size**) to be kept for downstream analysis.

Is it possible to use limma DE analysis also with groups composed by only 2 samples? NB. It is possible that in this particular dataset the smallest group size is 3.

If not, which alternative should I use?

### Update

I have 17 samples. My idea is to test one group vs the rest of groups (all together). So at least the reference group would be composed by \> 10 samples. In this case, what can be a feasible analysis for DE?

1. Smyth, G. K., Law, C. W., Alhamdoosh, M., Su, S. & Ritchie, M. E. RNA-seq analysis is easy as 1-2-3 with limma, Glimma and edgeR. F1000Research 5, 1408 (2016).

# Answer
> 3 votes
The problem with most of the methods is that use the gene's variance for each group, which can't be calculated (reliable) when the sample is \<= 2. Also, statistically, it would have extremely low power, so the conclusions couldn't be trusted much.

You could also calculate the "raw" fold change by yourself (ie without the variance estimation and adjusting). If you had one more sample per group you could use DESeq2 which was thought for such cases. But from the paper:

> However, if there are two or fewer replicates for a condition, these samples do not contribute to outlier detection, as there are insufficient replicates to determine outlier status.

Which makes me unsure if it will work well for just two samples as it won't help to detect outliers genes.

# Answer
> 3 votes
While I share Llopis' concern about estimating variance from 2 samples, the statement you quoted is about avoiding false positives from genes that are only expressed in a few samples. It's fairly common to exclude genes that expressed in fewer samples than the smallest group even if the number of samples per group is much higher than 3.

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


# Question
Title: How to split bedfile into non-overlapping regions and compute aggregation function on duplicate segments
I have a bed-like file in which I have several overlapping regions with associated features. I want to split the overlapping regions into disjoint regions in order to process the file (e.g. with pandas), group by regions and compute an aggregate function.

The input data format is the following (overlapping are marked with `*`):

```
chr9    8855925    8986425    -1  *
chr9    8855965    8986325     0  *
chr9    8893857    9022105    -1

```

I want to split into disjoint intervals, then I want to groupby per-interval and compute an aggregation function (e.g. the mean of the last field):

```
chr9    8855925    8855965    -1
chr9    8855965    8893857  -0.5
chr9    8893857    8986325 -0.67
chr9    8986325    8986425    -1
chr9    8986425    9022105    -1

```

I was thinking of using bedtools and pandas to aggregate. However, I was unable to find the right tool to do this job in the bedtools/pybedtools documentation.

# Answer
> 3 votes
To make disjoint intervals, you could use BEDOPS `bedops --partition`, piping to `bedmap --mean` to get the mean signal over disjoint regions.

Starting with the input bedgraph file, convert it to five-column BED with GNU `awk`, putting the signal in the fifth column per UCSC convention:

```
$ awk -vOFS="\t" '{ print $1, $2, $3, ".", $4 }' /tmp/in.bedgraph | sort-bed - > /tmp/in.bed
$ cat in.bed
chr9    8855925 8986425 .   -1
chr9    8855965 8986325 .   0
chr9    8893857 9022105 .   -1

```

The disjoint set from these intervals would look like this:

```
$ bedops --partition /tmp/in.bed
chr9    8855925 8855965
chr9    8855965 8893857
chr9    8893857 8986325
chr9    8986325 8986425
chr9    8986425 9022105

```

We can pipe these intervals to `bedmap --mean` with the original five-column BED file as a "map" file, calculating the mean signal over the disjoint intervals:

```
$ bedops --partition /tmp/in.bed | bedmap --echo --mean --delim '\t' - /tmp/in.bed > /tmp/answer.bedgraph
$ cat /tmp/answer.bedgraph
chr9    8855925 8855965 -1.000000
chr9    8855965 8893857 -0.500000
chr9    8893857 8986325 -0.666667
chr9    8986325 8986425 -1.000000
chr9    8986425 9022105 -1.000000

```

The `--mean` option calculates the arithmetic mean. Other aggregation functions are available to process the signal or score (fifth) column in map BED files, like `--stdev`, `--median`, `--mad`, etc.

When calculating statistics, the default precision is six digits. The `--prec <n>` option can be used to specify desired precision in the score result.

See `bedmap --help` or the online documentation for a full description.

To avoid making and cleaning up intermediate files, `bash` process substitutions could be used to generate a one-liner:

```
$ bedops --partition <(sort-bed in.bedgraph | awk -vOFS="\t" '{ print $1, $2, $3, ".", $4 }') | bedmap --echo --mean --delim '\t' - <(sort-bed in.bedgraph | awk -vOFS="\t" '{ print $1, $2, $3, ".", $4 }') > answer.bedgraph

```

This can be useful for cleaner pipelines or with processing very large, whole-genome scale inputs on systems with limited disk space.

---
Tags: bed, bedtools, pybedtools
---


# Question
Title: How to apply function on bedfile based on a reference bedfile
I have two bed-like files:

**File1**

```
chr1     1     10    .    1    Sample1
chr1    20    180    .    1    Sample1
chr1     1      5    .    1    Sample2
chr1    70    100    .    1    Sample3

```

**File2** (reference)

```
chr1     1     5     .    0
chr1     5    20     .    1
chr1    50   120     .   -1

```

I'd like to find a way to apply a function (in my case a subtraction) between the first file and the second, using the second as a reference, and producing a results file containing also the indication of the sample name. For example, in my case the results file would be:

```
chr1     1      5    .    1    Sample1
chr1     5     10    .    0    Sample1
chr1    20     50    .    1    Sample1
chr1    50    120    .    2    Sample1
chr1   120    180    .    1    Sample1
chr1     1      5    .    1    Sample2
chr1    70    100    .    2    Sample3

```

I was thinking about implementing it in pandas, but it would be long to do and probably not so efficient. I looked for bedmap but it seems to me that this operation can be carried on (with the sum operation) only within a single file during aggregation, not between two different bed files. So maybe it is not the right tool to use in this case.

Also, it would be important in the final file to retain the sample name.

# Answer
You could munge your first file into a BED5 file:

```
$ awk -vOFS="\t" '{ print $1, $2, $3, $6, $5 }' /tmp/f1.bed | sort-bed - > /tmp/f1.bed5

```

Sort the second file, if necessary:

```
$ sort-bed /tmp/f2.unknownSort.bed5 > /tmp/f2.bed5

```

Take the union of the two sorted files:

```
$ bedops -u /tmp/f1.bed5 /tmp/f2.bed5 > /tmp/union.bed

```

Here's what the union looks like:

```
chr1    1       5       .       0
chr1    1       5       Sample2 1
chr1    1       10      Sample1 1
chr1    5       20      .       1
chr1    20      180     Sample1 1
chr1    50      120     .       -1
chr1    70      100     Sample3 1

```

We have one name scheme for reference element IDs (`.`) and a second name scheme for sample IDs (`Sample*`).

We will exploit this difference in name schemes when we apply operations in a certain order. For this example, we will subtract a reference element score from a sample element score.

We can partition the union, mapping the disjoint elements back to the unioned set:

```
$ bedops --partition /tmp/union.bed | bedmap --echo --echo-map-id --echo-map-score --delim '\t' - /tmp/union.bed
chr1    1   5   .;Sample2;Sample1   0.000000;1.000000;1.000000
chr1    5   10  Sample1;.   1.000000;1.000000
chr1    10  20  .   1.000000
chr1    20  50  Sample1 1.000000
chr1    50  70  Sample1;.   1.000000;-1.000000
chr1    70  100 Sample1;.;Sample3   1.000000;-1.000000;1.000000
chr1    100 120 Sample1;.   1.000000;-1.000000
chr1    120 180 Sample1 1.000000

```

Each of the fourth and fifth columns are semi-colon-delimited strings. Each element is ordered by the order in which it is discovered in a `bedmap` operation.

Both columns have the same ordering or same one-to-one correspondence of ID to signal score. This is very useful, because we can turn these into a pair of ordered arrays and apply some logic, depending on the answer we want back.

For example, here is some `awk` logic to parse those two strings and subtract a reference element score from a sample element score:

```
$ bedops --partition /tmp/union.bed \
    | bedmap --echo --echo-map-id --echo-map-score --delim '\t' - /tmp/union.bed \
    | awk -vFS="\t" -vOFS="\t" \
        '{ \
             i=split($4,a,";"); \
             j=split($5,b,";"); \
             if ((i==1) && (a[i]==".")) { pass; } \
             else { \
                 for (k=1; k<=i; k++) { \
                     if (a[k]==".") { refIdx=k; } \
                 } \
                 for (k=1; k<=i; k++) { \
                     if (k!=refIdx) { print $0, ".", b[k]-b[refIdx], a[k]; } \
                 } \
             } \
         }' \
    | cut -f1-3,6-8 \
    | sort-bed - \
    > answer.bed

```

The file `answer.bed` will look like this:

```
chr1    1   5   .   1   Sample1
chr1    1   5   .   1   Sample2
chr1    5   10  .   0   Sample1
chr1    20  50  .   1   Sample1
chr1    50  70  .   2   Sample1
chr1    70  100 .   2   Sample1
chr1    70  100 .   2   Sample3
chr1    100 120 .   2   Sample1
chr1    120 180 .   1   Sample1

```

This partitioning of intervals is different from the one in your question. Perhaps this is due to how you partitioned your set. In any case, perhaps this would give you an idea of how to apply custom operations over a disjoint partitioning.

The main idea is that this part of the `awk` script finds the score associated with the reference file ("File2"):

```
for (k=1; k<=i; k++) {
    if (a[k]==".") { refIdx=k; }
}

```

Any index other than `refIdx` corresponds to a `Sample*` element. That means we can do operations between `Sample*` scores and the specific reference score.

We also leave out mappings where a reference score does not associate with any `Sample*` element:

```
if ((i==1) && (a[i]==".")) { pass; }

```

Thus, all the operations are always between a pairing of a sample and the mapped reference element.

The `cut` and `sort-bed` statements at the end just clean up the result so that you get something that looks like what you are expecting at the end, which is also sorted correctly and ready to use for downstream set operations.
> 2 votes

---
Tags: bed, bedtools
---


# Question
Title: Retrieving gene matrix from internet in GSEA desktop
This error appeared to me, and I have the program connected to the internet. I already installed java 8, but it still doesn't get me the gene set databases.

# Answer
Have you checked firewall settings on your laptop?

See **Firewall / FTP connection issues** at:

GSEA wiki - known issues
> 2 votes

---
Tags: rna-seq, gsea
---


# Question
Title: Raw vs Filtered in the output of cellranger count
After running `cellranger count` I got two relevant for further analysis folders: `filtered_gene_bc_matrices` and `raw_gene_bc_matrices`. What is the difference between them? What does `cellranger` filter out, and what should we preferably use for further analysis?

# Answer
I'm unsure whether this is the answer you are looking for, but when looking into 10X cellranger documentation for the Matrices Output:

> **Unfiltered gene-barcode matrices**: *Contains every barcode from fixed list of known-good barcode sequences. This includes background and non-cellular barcodes.*
>
> **Filtered gene-barcode matrices**: *Contains only detected cellular barcodes.*

Source
> 5 votes

---
Tags: scrnaseq, 10x-genomics
---


# Question
Title: How to extract the protein fasta file from a genbank file?
Let's say I have a genbank file, e.g. the one downloaded from here, that contains entries as

```
 gene            190..255
                 /gene="thrL"
                 /locus_tag="b0001"
                 /gene_synonym="ECK0001; JW4367"
                 /db_xref="EcoGene:EG11277"
                 /db_xref="GeneID:944742"
 CDS             190..255
                 /gene="thrL"
                 /locus_tag="b0001"
                 /gene_synonym="ECK0001; JW4367"
                 /function="leader; Amino acid biosynthesis: Threonine"
                 /GO_process="GO:0009088 - threonine biosynthetic process"
                 /codon_start=1
                 /transl_table=11
                 /product="thr operon leader peptide"
                 /protein_id="NP_414542.1"
                 /db_xref="ASAP:ABE-0000006"
                 /db_xref="UniProtKB/Swiss-Prot:P0AD86"
                 /db_xref="EcoGene:EG11277"
                 /db_xref="GeneID:944742"
                 /translation="MKRISTTITTTITITTGNGAG"

```

and I want to extract the potein fasta from it, how would I do this?

I found a code sample here. When I then run:

```
from Bio import SeqIO

gbk_filename = 'NC_000913.3.gb'
faa_filename = 'NC_000913.3_warwick.fasta'
input_handle  = open(gbk_filename, "r")
output_handle = open(faa_filename, "w")

for seq_record in SeqIO.parse(input_handle, "genbank") :
    print "Dealing with GenBank record %s" % seq_record.id
    for seq_feature in seq_record.features :
        if seq_feature.type=="CDS" :
            assert len(seq_feature.qualifiers['translation'])==1
            output_handle.write(">%s from %s\n%s\n" % (
                   seq_feature.qualifiers['locus_tag'][0],
                   seq_record.name,
                   seq_feature.qualifiers['translation'][0]))

output_handle.close()
input_handle.close()

```

it results in

```
     31     for seq_feature in seq_record.features :
     32         if seq_feature.type=="CDS" :
---> 33             assert len(seq_feature.qualifiers['translation'])==1
     34             output_handle.write(">%s from %s\n%s\n" % (
     35                    seq_feature.qualifiers['locus_tag'][0],

KeyError: 'translation'

```

How could this be fixed?

# Answer
One can get it to work by using `SeqIO.InsdcIO.GenBankCdsFeatureIterator`:

```
from Bio import SeqIO

file_name = 'NC_000913.3.gb'

# stores all the CDS entries
all_entries = []

with open(file_name, 'r') as GBFile:

    GBcds = SeqIO.InsdcIO.GenBankCdsFeatureIterator(GBFile)

    for cds in GBcds:
        if cds.seq is not None:
            cds.id = cds.name
            cds.description = ''
            all_entries.append(cds)

# write file
SeqIO.write(all_entries, '{}.fasta'.format(file_name[:-3]), 'fasta')

```

This will then give the desired output (just the first two entries):

```
>b0001
MKRISTTITTTITITTGNGAG
>b0002
MRVLKFGGTSVANAERFLRVADILESNARQGQVATVLSAPAKITNHLVAMIEKTISGQDA
LPNISDAERIFAELLTGLAAAQPGFPLAQLKTFVDQEFAQIKHVLHGISLLGQCPDSINA
ALICRGEKMSIAIMAGVLEARGHNVTVIDPVEKLLAVGHYLESTVDIAESTRRIAASRIP
ADHMVLMAGFTAGNEKGELVVLGRNGSDYSAAVLAACLRADCCEIWTDVDGVYTCDPRQV
PDARLLKSMSYQEAMELSYFGAKVLHPRTITPIAQFQIPCLIKNTGNPQAPGTLIGASRD
EDELPVKGISNLNNMAMFSVSGPGMKGMVGMAARVFAAMSRARISVVLITQSSSEYSISF
CVPQSDCVRAERAMQEEFYLELKEGLLEPLAVTERLAIISVVGDGMRTLRGISAKFFAAL
ARANINIVAIAQGSSERSISVVVNNDDATTGVRVTHQMLFNTDQVIEVFVIGVGGVGGAL
LEQLKRQQSWLKNKHIDLRVCGVANSKALLTNVHGLNLENWQEELAQAKEPFNLGRLIRL
VKEYHLLNPVIVDCTSSQAVADQYADFLREGFHVVTPNKKANTSSMDYYHQLRYAAEKSR
RKFLYDTNVGAGLPVIENLQNLLNAGDELMKFSGILSGSLSYIFGKLDEGMSFSEATTLA
REMGYTEPDPRDDLSGMDVARKLLILARETGRELELADIEIEPVLPAEFNAEGDVAAFMA
NLSQLDDLFAARVAKARDEGKVLRYVGNIDEDGVCRVKIAEVDGNDPLFKVKNGENALAF
YSHYYQPLPLVLRGYGAGNDVTAAGVFADLLRTLSWKLGV

```
> 6 votes

---
Tags: fasta, python, genbank
---


# Question
Title: Can I run STAR without an annotation file?
I wish to use Rascaf to scaffold a fragmented draft genome.

For this, I need to provide a BAM file of aligned RNA-seq reads and the draft genome.

So, I indexed the draft genome with STAR like this:

```
STAR --runMode genomeGenerate --genomeDir output/index --genomeFastaFiles draft.fa

```

Then, I tried to aligned the reads like this:

```
STAR --genomeDir output/index --readFilesIn reads.fastq.gz --readFilesCommand gunzip -c --outFileNamePrefix output/alignment --quantMode GeneCounts --outSAMunmapped None --outSAMtype BAM SortedByCoordinate --outSAMattrRGline ID:RG1 CN:yy \"DS: z z z\" SM:sample

```

The latter commands outputs the following error:

```
Transcriptome.cpp:51:Transcriptome: exiting because of *INPUT FILE* error: could not open input file exonGeTrInfo.tab
Solution: check that the file exists and you have read permission for this file
          SOLUTION: utilize --sjdbGTFfile /path/to/annotantions.gtf option at the genome generation step or mapping step

```

It suggests I should use a GTF file. But this is a draft genome of an individual for which no GTF is available. I could try to adapt a GTF from a close relative but for the scaffolding, it seems superfluous. Is there a way I can map with STAR without the GTF?

These commands were run on a GNU/Linux machine.

# Answer
I don't think you can use the `--quantMode GeneCounts` option with no annotations. I think the error is trying to look for an exon file generated from the annotations to do the quantitation on. Remove that and I think it should work, as the manual specifically states that annotations are optional but highly recommended.
> 3 votes

---
Tags: read-mapping, rna-alignment, genome-sequencing, scaffold
---


title: Why Ti/Tv ratio?

text:
# Question
Title: Why Ti/Tv ratio?
I'm interested in the transition/transversion (Ti/Tv) ratio:

> In substitution mutations, transitions are defined as the interchange of the purine-based A↔G or pryimidine-based C↔T. Transversions are defined as the interchange between two-ring purine nucleobases and one-ring pyrimidine bases.

(Wang *et al*, *Bioinformatics*, **31** (3), 2015)

* How exactly does this ratio imply false positives? Too high ==\> high false positive rates? Or too low?
* Why is the expected value for random substitutions for the Ti/Tiv ratio 0.5?
* If the ratio is expected to be 2.10 for WGS, but I get 3.00 what does that mean? What if I get 1.00?

# Answer
> 8 votes
I more like to use "ts/tv" for transition-to-transversion ratio. This abbreviation had been used in phylogenetics. When NGS came along, some important developers started to use "ti/tv", but I am still used to the old convention.

> Why is the expected value for random substitutions for the Ti/Tv ratio 0.5?

There are six types of base changes. Two of them are transitions: `A<->G` and `C<->T` and the other four types are transversions. If everything were random, you would expect to see twice as many transversions – ts:tv=2:4=0.5.

> If the ratio is expected to be 2.10 for WGS

The expected ratio is not "2.10 for WGS". It is 2–2.10 for human across the whole genome. You see this number when you align the human genome to chimpanzee or when you focus on an accurate subset of human variant calls. However, in other species, the expected ts/tv may be very different. Also, this number is correlated with GC content. You get a higher ts/tv in high-GC regions, or in coding regions which tend to have higher GC, too. Partly as a result, it is hard to say what is expected ts/tv in exact.

> but I get 3.00 what does that mean? What if I get 1.00?

If you get 3.00, your callset is highly biased. If you get 1.00, your callset has a high error rate. Suppose the exact ts/tv is $\beta$ and you observe ts/tv $\beta'\le\beta$, you can work out the fraction of wrong calls to be (assuming random errors have ts/tv=0.5) $$ \frac{3(\beta-\beta')}{(1+\beta')(2\beta-1)} $$ This is of course an approximate because $\beta$ is not accurate in the first place and because errors are often not random, so their ts/tv is not really 0.5.

> How exactly does this ratio imply false positives? Too high ==\> high false positive rates? Or too low?

Too low ==\> high false positive rate; too high ==\> bias. In practice, you rarely see "too high" ts/tv.

# Answer
> 3 votes
I haven't done variant calling myself, but 5 seconds of googling lead me to this:

> This metric is the ratio of transition (Ti) to transversion (Tv) SNPs. If the distribution of transition and transversion mutations were random (i.e. without any biological influence) we would expect a ratio of 0.5. This is simply due to the fact that there are twice as many possible transversion mutations than there are transitions. However, in the biological context, it is very common to see a methylated cytosine undergo deamination to become thymine. As this is a transition mutation, it has been shown to increase the expected random ratio from 0.5 to ~2.01. Furthermore, CpG islands, usually found in primer regions, have higher concentrations of methylcytosines. By including these regions, whole exome sequencing shows an even stronger lean towards transition mutations, with an expected ratio of 3.0-3.3. A significant deviation from the expected values could indicate artifactual variants causing bias. If your TiTv Ratio is too low, your callset likely has more false positives.

Since methylated cytosine deamination is a fairly universal phenomenon, if you're seeing a ratio that deviates substantially from this expectation, you've either discovered some exciting new biology, or there are issues with your variant calling pipeline.

---
Tags: variant-calling
---


# Question
Title: Difference between copy number neutral reads and active reads
In this paper, the authors talk about *copy number neutral reads* (as reads that equally represent parental chromosomes) and *active reads* (as reads from only one parent chromosome):

> We reasoned that the reads mapped into a genomic window can be partitioned in two sets: one set includes reads that equally represent parental chromosomes (copy number neutral reads); and the other set contains reads from only one parent chromosome (active reads).

Nonetheless, I have some difficulties understanding these two concepts. Can somebody try to expand on what they wrote?

# Answer
In genetics normally one gets a copy (allele) of a gene (or actually a chromosome) from each parent. If the alleles differ between parents, for example by a SNP, one would have one copy of each allele. With (NGS) sequencing you would get an even amount of reads for both copies. However with somatic mutations, hence mutations which occur after meiosis and gamete forming, one can lose or gain parts of the chromosome (copy number variation). In this situation, some genomic regions might be deleted or amplified. So the allele of one parent can be deleted, leaving you with only the allele of the other parent.

As requested, I have made a simple drawing of a situation where a region differs between parent alleles with SNPs. After a deletion of one of the allele only reads of the other allele are left.
> 1 votes

---
Tags: ngs, genomics, cancer
---


# Question
Title: Motif Discovery in Bipartite networks using Cytoscape
I am trying to find over-expressed motifs in a bipartite transcription factor gene network.My genes are not very well characterized, so finding functional motifs has been a challenge. I used the motif-discovery app on Cytoscape and calculated the total number of 3 node motifs on the network. In my network, where transcription factors regulate genes, the only kind of three node motif that is possible is a simple input module.

This is an example of a simple input module where a transcription factor regulates multiple genes

I further created random networks to my real world network and statistically (Student's t-test) proved that there are more number of single input modules in the real rather than the random network. However, I am just not feeling satisfied with just this one explanation because I would have like to draw some biological inference from this. I was thinking of also analysis the 4 node and 5 node motifs in the network, find the number of times any of my hub nodes occupy a single node of these motifs and maybe use that and prove the importance of my hub nodes (though hub nodes are worthy of more exploration just by virtue of their high degree). Has somebody done something like this previously. It is a given that hub nodes would be part of the maximum number of motifs just by virtue of their high degree and connection to multiple genes. So would this prove to be redundant?

Thank you in advance for any help that could be provided.

# Answer
I have recently develop an R package (sorry not a Cytoscape solution) to deal with this kind of situations of something grouping other elements, like pathways and genes.
The package is based on GSEABase of Bioconductor to store this relationships. Basically you should store a gene set for each TF and genes it is linked. I would recommend to use the name of the TF as the name of the gene set and then create a GeneSetCollection to be able to use the package.

The package provides with functionalities that allow you to calculate the distribution of genes per gene set (using `genesPerPathway`), as well as, `pathwaysPerGene` to know in how many other TF regulate each gene.

I am not sure of the biological inference you can take from this (I designed the package to test this) but you can also see if there is some specific preference for the TF to link to genes which are also regulated by other TF. This could be done by looking at the conditional probability (`condPerPathways` or `sizesPerPathway`).
Sorry the manual pages don't have examples and the vignettes are incomplete yet.

Also you could use the package tidygraph to calculate the centrality of the TF or of the genes you want.
> 2 votes

---
Tags: networks, motifs
---


# Question
Title: How do I find a list of PDB structures bound to a peptide inhibitor or peptide-like inhibitor?
I looked up existing databases but ones like iPPI-DB or enter link description here either provide non-peptide inhibitors or only sequences. I need structures of proteins bound to other proteins, the latter of which is an inhibitor of the former, and the complex has a structure deposited in the PDB database. Here's an example: PDB:3U1J

In case there is no existing database in literature, it'd be great even if I can screen for such PDB entries using some code, but I can't come up with an algorithm for it. Thank you!

# Answer
The quickest way to screen the PDB database is *Advanced Search* from either RCSB or PDBe.

Since you need entries with 2+ proteins you can set the number of protein chains to at least 2 in the RCSB search. Or you can select "protein/protein complex" in *Assembly composition* in the PDBe search.

It's less clear how to search for inhibitors. If you'd look for the "inhibitor" word in the molecule names you'd get (together with the previous criterium) 1200+ entries.
> 1 votes

---
Tags: pdb
---


# Question
Title: How to download SNP data from specific regions, population and positions?
I am totally new in Bioinformatics and I would like to apply my knowledge in feature selection on the tag SNP problems.

To do that, I've read a lot of papers and books in order to understand the main concepts and how the algorithms to this task work.

However, is been very hard for me to understand the organization of the huge amount of genomic data available over the internet.. HapMap, 1000 Genomes, dbSNP and so on... (!) I downloaded GB of data from different sources, but I really do not know what to do with then.

--

To be more specific in my question:

First I would like to run a Genetic Algorithm to find the tag SNP based on a simple linkage disequilibrium coverage measure. The authors of this GA used the following data

> Three published SNP data sets were downloaded from HapMap (http://www.hapmap.org/) ***(NOT AVAILABLE ANYMORE)*** for evaluating the prediction accuracies achievable by our method:
>
> • ENr112. 343 SNPs genotyped in population ASW (African ancestry in Southwest USA) on chromosome 2 from position 51512208 to 52012208. The genotypes of 83 individuals were used.
>
> • ENr113. 380 SNPs genotyped in population ASW on chromosome 4 from position 118466103 to 118966103. The genotypes of 83 individuals were used.
>
> • ENm013. 683 SNPs genotyped in population ASW on chromosome 7 from position 89621624 to 90736048. The genotypes of 83 individuals were used.

Besides that, the authors stated that

> We represent a genotype g\_i of length L by a sequence over {0, 1, 2}^L (since we assume only bi-allelic SNPs). 0 and 1 stand for the minor and major homozygous types respectively, meaning that the two alleles at that SNP locus are either 0 or 1; the major and minor alleles are determined by the frequency at which each allele appears in the population. A value of 2 denotes the heterozygous type for an SNP having two different alleles.

And also

> To obtain phasing of genotypes, we used the GEVALT algorithm.

I believe they obtain the aforementioned data in genotype format (something with 0, 1 and 2 values only) and then phasing it to apply in the problem. The image below is a example of this process.

So

1. How and where can I download these data?
2. I know HapMap is not available due to a security issue but there are the FTP, but the FTP folders are really complex to understand. Can someone guide me to download data from there?

I hope I made myself clear.

# Answer
The HapMap website was suspended due to a security issue. Archived HapMap data is available via FTP, but it's a better idea to download up-to-date International Genome data from the IGSR website (which includes the HapMap individuals as a subset).

The format you use will be dependent on the data analysis program you're using. In most cases involving SNP analysis, the VCF format is probably what would be used.
> 2 votes

---
Tags: snp, public-databases, genomics, data-download, haplotypes
---


# Question
Title: Network comparison of single cells (from sequencing data)
What methods can we use to compare networks (PPI, gene regulatory etc) created from single-cell gene expression (or proteomics) data?

There are methods that construct networks by comparing two conditions (e.g. from (gene) co-expression, or differential gene correlation (DGCA, Differential Gene Correlation Analysis)), and infer properties based on these networks.

What I'm interested here is the comparison of networks constructed for each single cell. The networks can be constructed solely from the data, or using known interactions (edges) with the help of a database. These comparisons could be used for understanding or classifying cells.

This will probably require a network construction method (which may be hard because we do not have replicates for single cells), some form of network similarity measure / network alignment method. This is a very hard task in general, but here we have 1-to-1 correspondence for each node between the cells (networks), which should help.

Similar discussions at

One related paper for single-cells is Chan et al., 2017, which uses groups of cells (created by other methods, e.g. clustering) to construct networks, and manually compares them (common edges / Venn diagram).

# Answer
You can't build a network of a single cell only with the expression of a single cell. You either need previous known interactions or pathways or you need to use several cells/samples. If you use previous known information, you can use pathway information, otherwise you can group some cells and use something along the lines of WGCNA to find a scale-free network.

In the first case one can add the expression upon it and see what is happening, what edges are kept and which ones not...

---

Once you have the network for each cell you can compare the metrics for each cell. Some metrics to look at are the centrality of each gene and the betweneess.
In a more general way, you can look if it follows a scale-free topology or other type of network structures, or how many nodes are present, how many edges are present...
> 4 votes

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


# Question
Title: Expand granges object different amounts upstream vs. downstream
I am attempting to get gene regions and their immediate neighborhoods using Bioconductor GenomicRanges packages.

It is very easy for me to extend the regions around genes if I want to extend it the same distance upstream and downstream as in this example:

```
library(GenomicRanges)
library(Homo.sapiens)

genes = genes(Homo.sapiens)
genes_5kb_flank = genes + 5000

```

And if I want to extend both upstream and downstream different amounts centered on the TSS I can use `promoters(genes, upstream = 5000, downstream= 1000)`

But I am having trouble finding a function that is strand-aware and also chromosome length aware that allows me to extend the gene region a different amount upstream vs. downstream... For instance if I wanted a 5kb upstream flanking region and a 1kb downstream flanking region as is used for genomic region annotations on GREAT.

Anyone know of a drop in function for this? something like `promoters()` but for the full gene region instead of only around the TSS?

# Answer
> 2 votes
There doesn't appear to be a built in function for that, so here's a modified version of `promoters()`:

```
expandRange = function(x, upstream=2000, downstream=1000) {
  strand_is_minus = strand(x) == "-"
  on_plus = which(!strand_is_minus)
  on_minus = which(strand_is_minus)
  start(x)[on_plus] = start(x)[on_plus] - upstream
  start(x)[on_minus] = start(x)[on_minus] - downstream
  end(x)[on_plus] = end(x)[on_plus] + downstream
  end(x)[on_minus] = end(x)[on_minus] + upstream
  x
}

```

You can then call `expandedGenes = expandRange(genes, 5000, 1000)`.

---
Tags: bioconductor, sequence-annotation, genomicranges
---


# Question
Title: RunUMAP in Seurat not working: module 'umap' has no attribute 'UMAP'
I am trying to run `UMAP` in the following way:

```
RunUMAP(seurat_object, dims.use = 1:15, reduction.use = "pca", reduction.name = "umap",
    reduction.key = "UMAP", n_neighbors = 30L, min_dist = 0.3, metric = "correlation")

```

Or just this:

```
RunUMAP(seurat_object, dims.use = 1:15, reduction.use = "pca")

```

It is giving me the error:

> Error in py\_get\_attr\_impl(x, name, silent) : AttributeError: module 'umap' has no attribute 'UMAP'

I already installed `umap-learn` successfully in this way: `conda install -c conda-forge umap-learn` and checked that `Rstudio` is using the same `python` version as is provided by default in the terminal: `3.6.5`

# Answer
> 1 votes
This is resolved by installing the more recent version of the code from GitHub.

```
git clone https://github.com/lmcinnes/umap
cd umap
pip install --user -r requirements.txt
python setup.py install --user

```

As Peter said, this kind of question is best directed to the Issues page on their GitHub repository.

---
Tags: scrnaseq, seurat, umap
---


# Question
Title: How to use variables while calling pybedtools
I am trying to use bedtools and fetch sequences from a whole genome fasta file inside the script I get region coordinates as variables. For example:

```
chr_ = "chr1"
start = 3000
end = 3402

```

I am not sure how to wrap my function, because the example provided by BedTool.sequence only uses coordinates in a string form.

Can you please help me to have something like this?

```
a = pybedtools.BedTool(chr_, start, end, from_variables=True)

```

# Answer
You have to take your variables and make them into a single string variable.

for example:

```
chr_ = "chr1"
start = 3000
end = 3402
as_str = ' '.join([chr_, str(start), str(end)])
# or
as_str2 = "{} {} {}".format(chr_, start, end)

```

Once you have the string built you can use the example in the bedTool docs.

```
a = BedTool(as_str, from_string=True)

```
> 2 votes

---
Tags: bedtools, pybedtools
---


# Question
Title: determining doublets in single-cell RNA-seq
Doublets are a known problem with scRNA-seq experiments, where 2 or more cells are sometimes captured instead. To determine their presence, there are studies that mix multiple species (such as human and mouse) or distinct cell types and then count how many cells contain genes from the very different transcriptomes. In those experiments, you can pick any cell and determine if it's a doublet based on a few key genes. In a typical single-cell experiment on one sample, you don't have the luxury of knowing what cells to expect and if a particular cell looks abnormal.

If you have a large fraction of doublets, you should be able to see a double peak if you plot the number of genes/transcripts per cell. However, if you have something like 1% doublets, you probably wouldn't notice anything. It seems that most people simply filter out top X% cells based on the number of genes/transcripts. How do you know what X should be? Can you actually estimate the fraction of doublets in a given experiment?

Update:

I found a nice figure from Stoeckius et al. to illustrate the issue:
Although the number of transcripts/UMIs seems to be roughly double for doublets/multiplets (median of ~1000 vs ~500), there is still high overlap.

# Answer
# Expected rates of doublets / duplets / multiplets

Rates in the *'Chromium Single Cell 3’ Reagent Kits v2 User Guide':*

```
Multiplet Rate (%)  # of Cells Loaded   # of Cells Recovered
0.40%   870 500
0.80%   1700    1000
1.60%   3500    2000
2.30%   5300    3000
3.10%   7000    4000
3.90%   8700    5000
4.60%   10500   6000
5.40%   12200   7000
6.10%   14000   8000
6.90%   15700   9000
7.60%   17400   10000

```

---

# Detection and identification of doublets and multiplets

## Experimental

* **Visual inspection** (e.g. Fluidigm C1)
* **FACS index sorting data** can be used to identify doublets.
* **Mixed species experiment**.
* **Cell "hashing" with oligo-tagged antibodies** (Stoeckius et al., 2017): cells (barcodes) positive for two or more oligos are multiplets.

## Computational

* Mixed gene expression (e.g. lymphoid and myeloid, or T and B cell)
* Ilicic et al. describe an approach for detecting low quality cells (incl. multiplets). They train a support vector machine (SVM) model on a dataset annotated by microscopy inspection. They find, among others, that

+ *"wells containing multiple cells (multiples) show similar expression and noise patterns to broken cells"*, suggesting *"that multiple cells contain a mixture of broken and high quality cells."*
    + difference in proportion of duplicated reads
    + higher ERCC spike-in ratios
* demuxlet (Kang et al., 2017, Nature Biotechnology)
* DoubletDetection
> 10 votes

# Answer
**We cannot assume that doublets will produce more UMIs**

I would caution the assumption that all doublets will have twice the UMI levels of isolated single-cells. Many "doublets" could contain multiplets of 3 or more cells, depending on how many cells have been loaded into the experiment. Most vendors of single-cell technologies quote multiplet rates for doublets, triplets, etc combined.

Furthermore, we cannot assume that the data from multiplet (doublet) cells is of high quality or produces comparable reads and UMI counts than we would expect from two single cells combined. Some tests of doublets (such as using mixed mouse and human samples) have found that the data generated from doublet droplets or nanowells has relatively low UMI abundance. This may be due to the protocols having been optimised for droplets or nanowells containing only one cell. Especially where tissue preparation is required this is a concern and is consistent with observations that doublet cells appear to be broken or contaminated with debris when they can be identified. This is highly dependent on the tissue you and working with, sample preparation, and the platform used to separate single-cells.

I have not been able to reproduce the above result using data from various species (using total UMI counts per droplet or housekeeper genes), although this is work still in progress. Even in the data shown, many multiplet samples have similar total UMI counts to the true single-cells. Therefore, they cannot be separated by a simple threshold (despite being significantly different). There have been several tools developed already for detecting multiplets (and filtering them out) but this is still an area of active development and there are more ways that these could be refined. In many single-cell analyses, arbitrary thresholds are being used to filter the data before further analysis and this may be inadvertently reducing the impact of multiplets if they have produced fewer UMI counts than single-cells for which the protocols have been optimised for.

**Experimental protocols can be optimised to reduce doublets**

If you have any insight into how these bioinformatics solutions to issues with multiplets could be used or improved I encourage you to pursue it. However, I suspect that the eventual solution to this issue is in optimising the single-cell NGS experiments for the sample type and platform that you are working with. Technology in single-cell 'omics is continuing to develop and any solution which produces more useful single-cell data than contaminants or doublets (that need to be removed), the more cells can be sampled at the same cost. There is a strong incentive for experimentalists and single-cell technology vendors to optimise this issue.
> 4 votes

---
Tags: rna-seq, scrnaseq, single-cell, quality-control, 10x-genomics
---


title: HIVdb program and FASTA sequence input

text:
# Question
Title: HIVdb program and FASTA sequence input
I am trying to use HIVdb(https://hivdb.stanford.edu/hivdb/by-sequences/) program. I have a HIV protease which I want to test with most commonly used PIs. But when I entered it in fasta sequence format to the HIVdb, it give me error code saying

> Critical: There were no Protease, Reverse Transcriptase, or Integrase genes found, refuse to process.

For example, if I download FASTA sequence of 1HPV from PDB and paste it in HIVdb input sequence window as below,

```
PQITLWQRPLVTIKIGGQLKEALLDTGADDTVLEEMSLPGRWKPKMIGGIGGFIKVRQYDQILIEICGHKAIGTVLVGPT
PVNIIGRNLLTQIGCTLNF

```

I get the error code. Any suggestion on how to input sequences in HIVdb? (Note: I am new to this subject)

# Answer
The HIVdb page states the following in the second paragraph (emphasis mine):

> Sequences can be entered as plain text if just one sequence is entered. Sequences must be entered using the FASTA format if multiple sequences are entered. Sequences can be pasted in the text box or uploaded using the File Upload option. The upper limit is currently 1000 sequences containing 3000 **nucleotides** per sequence.

The sequence in your question is a protein, so it consists of amino acids, not nucleotides. So, if you instead give a nucleotide sequence, for example this one:

```
>KX692212.1 HIV-1 isolate CH47 from Switzerland pol protein (pol) and gag protein (gag) genes, partial cds
CCTCAGATCACTCTTTGGCAACGACCCCTCGTCACAATAAAGATAGGGGGGCAACTAAAGGAAGCTCTAT
TAGATACAGGAGCAGATGATACAGTATTAGAAGAAATGAGTTTGCCAGGAAGATGGAAACCAAAAATGAT
AGGGGGAATTGGAGGTTTTATCAAAGTAAGACAGTATGATCAGATACTCATAGAAATCTGTGGACATAAA
GCTATAGGTACAGTATTAGTAGGACCTACACCTGTCAACATAATTGGAAGAAATCTGTTGACTCAGATTG
GTTGCACTTTAAATTTT

```

The program will work as expected:

---

For future reference, I found the KX692212 sequence above by taking your protein sequence and running a tBLASTn analysis against HIV sequences in NCBI's `nr` database. That was the first hit.
> 1 votes

---
Tags: fasta
---


# Question
Title: How to search a specific sequence in BAM files for 10X experiment
I have to search a specific sequence in a set of cells (\> 1k cells) from a single-cell experiment done with 10X genomics.

As input file I have a single bam file, and 24 fastqs, therefore each file contains information about several cells. How can I search for a specific sequence and underpin which single cells have it?

I should be able to search for a specific sequence in a bam using samtools, etc. But the issue here is how to demultiplex the data from different cells to understand which cell provide the sequence.

# Answer
> 2 votes
The bam file has bam tags which say what reads belong to what cells.

https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/output/bam

---
Tags: rna-seq, scrnaseq, 10x-genomics
---


# Question
Title: Generating random BED intervals given constraints
Problem is to generate a random BED interval given the following constraints:

1. minimum `start`
2. maximum `end`
3. fixed length
4. maximum number of masked bases (similar to `-maxN` option in `faSplit`)
5. set of intervals to avoid overlap with
6. stay within chromosome size constraints

I can think of computationally intensive ways to code it up from scratch, but am wondering if there is a more efficient approach before re-inventing the wheel.

```
seq_records = {x.name: x for x in SeqIO.parse('path/to/genome.fa', 'fasta')}

def generate_random_interval(chrom, lower, upper, length, maxrep=1.0, avoid_intervals=None):
    while True:
        start = np.random.randint(lower, upper - length)
        end = start + length
        regenerate = False

        for interval in avoid_intervals:
            if (interval.start < end) or (interval.end > start):
                regenerate = True
        if 1. * seq_records[chrom].seq[lower, upper].count('N') / length > maxrep:
            regenerate = True

        if not regenerate:
            break

    return chrom, start, end

```

# Answer
You can perhaps collapse requirements for (1), (2), and (6), so long as bounds from (1) and (2) fall within the bounds of (6).

To help with this, you could use `fetchChromSizes` from the Kent UCSC utilities to quickly get chromosomes and chromosomal bounds for your genome of interest.

Given a chromosome and a `maximumEnd` that falls within its bounds, uniformly sample from `[0, maximumEnd-minimumStart-maxLength)`, adding `minimumStart` to determine start position and `maxLength` as length to get the stop position.

To obtain the desired number of samples in an efficient way, you could pipe these "candidate elements" into `bedmap --echo-map-size` operations on repeatmasked regions (4) and `bedops -n 1` on "regions-to-avoid" (5) operations to do rejection sampling.

Perhaps something like:

```
$ REPEATMASK_THRESHOLD=123
$ SAMPLES=1000
$ someScriptThatGeneratesCandidateIntervalsWithinBounds \
    | sort-bed - \
    | bedmap --echo-map-size --echo --delim '\t' - repeatmaskedRegions.bed \
    | awk -vthreshold=${REPEATMASK_THRESHOLD} '($1<threshold)' \
    | cut -f2- \
    | bedops -n 1 - regionsToAvoid.bed \
    | head -${SAMPLES} \
    > qualifyingCandidates.bed

```

This will be relatively efficient, because once the `head` process at the end writes the number of `SAMPLES` to `qualifyingCandidates.bed`, it will send a `SIGINT` signal up the pipeline that will terminate the upstream processes. So you're only doing sampling to the extent needed.

Hope this helps!
> 2 votes

---
Tags: python, bed, bedtools
---


# Question
Title: What kind of analysis can be done with differential expression of transcription factors?
I have two different stem cell types and their respective gene expressions from RNAseq.

I noticed that there is differential expression in some of the genes coding for histones as well as transcription factors?

What kind of interesting downstream analysis might be available?

Is it feasible to survey ENCODE for downstream TF binding sites? Can anything interesting be investigated with regard to the histone gene data?

Thanks in advance.

# Answer
> 1 votes
Expanding my comment.

TF usually bind in the upper region of the genes it regulates, so you could see which genes are in the region below the binding site. Perhaps there are some databases for TF. I found a website base on the ENCODE project with the genes regulated by each TF.You can use them for pathway or gene set enrichment analysis, but I am not aware of any specific analysis for genes regulated by a single TF.

The histone data is harder if it was methylation I would say to look if different states of the histones are in different spatial state, but with gene expression is harder to conclude anything. Nevertheless you can check methylation studies to see if gene expression is associated with methylation changes.

# Answer
> 1 votes
It sounds like you’ve already done differential expression analysis (with genes or transcripts) and a pathway enrichment analysis to find these. Differential expression of these functional groups is an interesting finding in its own right but there are several further investigations that could follow, depending on your biological research interests.

Clustering the expression patterns of the differentially expressed genes would be able to show if there are relationships between them (such as inter-gene correlations) or relationships to other phenotype data. This is commonly performed with unsupervised hierarchical clustering and displayed by a heatmap. You could also look into whether these differentially expressed genes are expressed in a tissue-specific manner. There is plenty of reference material available from ENCODE, FANTOM4/5, or GTEx that could be used for comparison (as well as cancer datasets) and the Human Cell Atlas Project (currently underway).

Another approach is pathway structure. There are many hierarchical annotation of gene functions such as gene ontology (GO terms) that overlap and are ranked: identifying smaller more specific pathways for future investigations. No pathway annotation is perfect so it’s best to replicate important findings across pathway databases (ideally different expression datasets as well). These functional groups of genes aren’t unordered: the structural relationships between pathways are well characterised in some cases. This includes many gene regulatory networks involving transcription factors. These graph structures are supported on databases such as WikiPathways and Reactome. Network analysis is relatively unexplored compared to differential expression.

While this is a bioinformatics question, there’s always a need to consider experimental validation as a future direction. If you already have a hypothesis based on your analysis, it may be more appropriate to design an experiment or contact experimental collaborators. If you have access to further genomics data, you could also explore how these differences are regulated: CAGE-Seq for enhancers and lncRNA, ATAC-Seq, and 3C/4C for chromatin structure, or Bisulfite-Seq and ChIP-Seq for epigenetic marks. An integrative analysis of expression and epigenomics data would be especially compelling and would likely account for different activity of transcription factor binding sites. Unfortunately, this can’t be done without more ‘omics data available but it would greatly benefit investigations of histones if it is feasible.

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


# Question
Title: Expression of a gene in different groups
I would like to check the expression of a gene in different groups like Disease vs Normal samples. I want to make a plot out of that to check whether it is significant or not.

From this paper lncRNA I see that Figure 1C they used RPKM value. But, I'm using read counts data. My question is can I use log cpm to make a plot like that or should I need to use only RPKM?

# Answer
> 2 votes
You should *never* use RPKM. It’s simply obsolete in the age of paired-end sequencing, and has been replaced by FPKM (which is, strictly speaking, a synonym).

The linked blog post explains more generally the problems that measures such as FPKM and CPM suffer from. A more robust measure is the TPM (transcripts per million), which scales CPM by the lengths of the transcripts.

Lastly, none of these measures perform very well for comparing gene expression across replicates. Since you have read count data, you should run DESeq2 or edgeR to perform differential expression analysis across all genes, and go from there.

As the comment noted, the paper you cite does not represent scientific best practice. In fact, a competent reviewer would have rejected this analysis.

# Answer
> 1 votes
For visualization purposes, using log(cpm) is fine. But plot don't check if differences are significant or not, statistical tests do. You can certainly add the results of a proper statistical test to a plot - like adding asterisks to denote a significant difference, for example.

# Answer
> 0 votes
You can use the raw counts to normalise your data. It is import to do so before doing a differential expression analysis: otherwise you will detect more reads from samples with higher coverage and from longer genes. A gene with longer mRNA molecules will produce more reads than a shorter one for the same number of transcript molecules. Samples will not produce the same number of reads due to variation in cDNA extraction, amplification, library preparation, etc.

This does not inherently mean that a sample expresses every gene more than another sample. Hence we compare samples by relative rather than absolute expression (for microarray and RNA-Seq data). This is in principle, very similar to using housekeeper genes in qPCR experiments except we adjust for total RNA or reads, gene length, and batch effects.

There are many ways to achieve this but I recommend using the *limma* package in R. This provides the “voom” normalisation procedure. “voom” is documented in a separate follow-up paper to the limma paper. It produces normalised data (on a log-scale) which is compatible with downstream analysis provided by the package (including built-in Venn diagram and volcano plot functions). See the package vignette for more details.

Of course, everything here only applies to bulk RNASeq experiments and procedures for single-cell analysis are different.

---
Tags: rna-seq, gene, differential-expression, rna, fpkm
---


# Question
Title: How to identify gene expression signatures from gene expression data?
I have TCGA gene expression data. I'm interested in identifying gene expression signatures using the data.

I would like to know whether there are any tools or R packages for identifying gene signatures.

How are gene signatures different from GSEA?

# Answer
> 7 votes
There is no golden/standard way to define a gene signature but they are completely different from a gene set enrichment analysis (GSEA). I will start with how to obtain a **gene signature**:

Usually this kind of signatures are defined by comparing a group against another ie doing a differential gene expression analysis. One then selects the most up or down regulated genes and defines them as the signature of the group of cells. That's how the gene sets of the MSigDB which are not from (or derived from) pathway databases or gene ontologies are defined. There is no gold threshold of a fold change or a specific number of genes, it (highly) depends on the experiment you run and how good are the isolation of the cells.
Tools: DESeq2, limma, edgeR,...

As suggested by gringer, another way is to find which genes are expressed together (independently of how do they relate to other groups). This method relays on having a good enough number of samples from the same condition to be able to identify which genes are expressed the same way in those cells. If these genes are not expressed that way in other cell lines, then you have your signature for this group. Again, there isn't a golden method, different types of co-expression networks can be build and depends on how specific do you want to be or which is the data provided.
Tools: WGCNA, diffcoexp, TReNA,...

This is different from a **gene *set* enrichment analysis**1, which consists on evaluating if a group of genes (usually a gene set, a pathway or genes from a specific gene ontology) are not randomly sorted in a sorted list of genes (usually this sorted list is the fold-change of a differential analysis).

There are several methods developed depending on the null hypothesis, the type of sub-method and what do they take into account...
Tools: GSEA, PLAGUE, GSVA, limma, fgsea,...

TLDR: To perform a gene set enrichment you don't need a gene signature but a gene set (that *can* be a gene signature or a pathway or a random group).

1: Note that gene set enrichment analyisis is different from gene enrichment analysis. For a good review on the methods I recommend this article

# Answer
> 4 votes
If you're interested in an unbiased approach you could have a look at Gene Correlation Network Analysis. WGCNA is a fairly old R package, with a somewhat confusing workflow and results that can be difficult to interpret. It also needs a large sample set, preferably from fairly different biological conditions, to produce meaningful results.

Asking how "gene expression signatures" differ from GSEA is like asking how surface textures differ from trademarks. Gene set enrichment analysis is for looking at gene rankings and seeing how they compare to a known / existing gene set. A gene set is considered "enriched" if the majority of the gene set appears in high-ranked genes (or alternatively in low-ranked genes), rather than a more even distribution of the set across the ranked gene list. Most commonly people will rank genes based on expression (or fold change), but the method is agnostic to the nature of the ranking.

# Answer
> 2 votes
Gene Signatures are sets of genes to report a consensus signal or function. They are not defined in a particular way but can be defined by any chosen group depending on the biological function that you are interested in. Typically these are defined by differential gene expression analysis or principal component analysis, although other approaches such as network analysis or time-course data could also be used.

There are plenty of examples of sets of genes published as gene signatures (especially in the field of cancer research). Those that have been cross-validated in other datasets or replicated by other research groups are ideal. If you wish to generate such a “signature” for a particular biological function, you could follow the approaches demonstrated in these published signatures (especially those that have been widely used and cited).

A gene expression signature is used to generate a consensus expression signal for these groups of genes (representing the given biological function). A widely accepted way to do this across samples of bulk RNASeq (or microarray data) is a “meta-gene” from the first eigenvector (v\[,1\]) of the singular value decomposition (SVD). This dimension reduction identifies a major contributor to variation across samples (in these genes) without negatively correlated genes averaging out. Dimension reduction can also be achieved by PCA and tSNE, depending on the requirements of your analysis.

---
Tags: rna-seq, r, gene, bioconductor
---


# Question
Title: Statistical approach to compare the SNP genotyping data among set of individuals
So, I have the genotyping data of about 650,000 SNPs for 96 individuals. I already know the Y DNA haplogroup of these individuals, so to some extent, I have a gross understanding of their ancestry.

What would be the best way to go about doing this? From what I have been reading, maybe a PCA would be a good idea? Suppose individuals 1-10 belong to a particular population (unique ancestary) then their SNP genotypes should be significantly similar to each other than say individuals 11-20 who belong to another population group. So, in theory I should get two clusters in my PCA plot.

I have the genotype of each SNP for each individual. How would I go about creating a PCA plot or any other alternative way of analysis for this kind of questions.

# Answer
PCA sounds like a good start to see a rough population structure among your 96 individuals, this step is half a sane thing to do, but also it can give you a hint how many distinct populations there are in your dataset. I have not personally run PCA on SNP data, but it seems that there are quite lot of tutorials online. For instance this one looks nice.

Once you know how many categories you can expect in your population, you can use Structure to classify individuals to the categories.
> 3 votes

---
Tags: snp, phylogenetics, snp-chip
---


# Question
Title: Is there any tools that can generating FBA model from genome sequences or proteome sequences?
Suppose I have all sequences of genes or proteins of certain species, are there any free tools I can use to generate a flux balance analysis model out of them?

What I was thinking is more like to annotate all the enzymes based on some annotation method. Then mapping these enzymes to known meta-metabolic networks that we know for instance the meta-metabolic network in KEGG database.

# Answer
> 2 votes
To generate flux balance, one should now the stoichiometric numbers of all the reactions that are taking place, as well as the conditions (that could affect the reactions). I am not sure we have enough knowledge, time and computationally power to do this.

At the moment I'm only aware of a paper simulating the *Mycoplasma genitalium* cell. But I'm not sure it was worth the effort as I haven't seen replicated or redone with other species, perhaps the model was not accurate enough to be worth it.

The tools used are here. The licenses are not a standard, but at least they seem free.

---
Tags: sequence-analysis, sequence-annotation, modelling
---


title: Viral Metagenomics

text:
# Question
Title: Viral Metagenomics
I am analyzing **viral metagenomics** data (Illumina Miseq) for the first time. I have used **Ray** (reference below) for *de novo* viral genome assembly before but I haven't done metagenomics analysis before.

I know that there are some tools like Metavelvet for metagenomics. But I still need to find time to learn those tools.

For now, can I do quick and dirty analysis using following method ?

1. De novo assemby using Ray
2. Create a personal database of all the viral ref genomes using NCBI database
3. Blast all the contigs generated by Ray against personal database.

I assume this should give me some idea about the viral genomes present in my reads. What could be the potential flaws of this method in comparison to the tools like Metavelvet that are dedicated to metagenomics?

---

Ray Meta: scalable de novo metagenome assembly and profiling. Sébastien Boisvert, Frédéric Raymond, Élénie Godzaridis, François Laviolette and Jacques Corbeil. Genome Biology (BioMed Central Ltd). 13:R122, Published: 22 December 2012

# Answer
The specifics will depend on the experimental setup and data, but a few general comments...

1) Metagenomic data is often large volume and high in microbial diversity but low and uneven in coverage; single genome assemblers may suffer particularly since abundant (or closely related) species/strains will tend to be interpreted as repeats within a genome rather than deriving from different genomes, leading to assembly errors. Checkout CAMI for some recent benchmarks of mategenomic assmeblers.

2) You might want to think about filtering your database to certain viruses, e.g. you don't expect to find RNA viruses in DNA sequencing. If the sample is an animal body site (like human gut) then consider restricting to viruses found in that organism (if it's an environmental sample this may not be possible).

3) BLAST will be slow, consider DIAMOND for a drop in faster replacement for BLASTP. Even faster approaches use only specific marker genes (like Metaphlan2 or Mocat2) or k-mer methods (like kraken or clark). Again CAMI benchmarks many of these and other metagenomic taxonomic classifiers.

Otherwise, I think possibly a missing step is filtering the non-viral sequences. You haven't said if your sequencing was enriched for viruses or depleted for bacteria/host prior to sequencing but if it isn't then the vast majority of your sequences will likely be non-viral (e.g. bacterial or host). Removing or marking these early can speed up downstream analysis and remove false positives (e.g. the human genome contains integrated viruses, bacteria contain prophages etc.). Also it is normally important to look for and remove potential contaminants like lab reagents or cloning vectors.
> 2 votes

---
Tags: assembly, covid-19, metagenome, illumina
---


title: Error creating indices using STAR

text:
# Question
Title: Error creating indices using STAR
I am trying to index wheat genome using STAR through following command

```
STAR --runMode genomeGenerate --genomeFastaFiles Triticum_aestivum.TGACv1.dna_sm.toplevel.fa --runThreadN 28

```

But getting following error,

```
terminate called after throwing an instance of 'std::bad_alloc'
what():  std::bad_alloc
Aborted (core dumped)

```

The size of the genome is 13.5 GB and the server capacity is 32 cores with 125 GB RAM, how this error can be resolved?

Memory usage:

I have decreased the number of threads to 8 but again getting the same error as shown below:

```
> STAR   --runMode genomeGenerate   --runThreadN 8   --genomeDir wheat   --genomeFastaFiles wheat.fa
Jun 04 02:13:31 ..... started STAR run
Jun 04 02:13:31 ... starting to generate Genome files
terminate called after throwing an instance of 'std::bad_alloc'
what():  std::bad_alloc
Aborted (core dumped)

```

Even I have tried to use STAR without any threads but still the same error:

```
> STAR   --runMode genomeGenerate   --genomeDir wheat   --genomeFastaFiles wheat.fa
Jun 04 02:19:44 ..... started STAR run
Jun 04 02:19:44 ... starting to generate Genome files
terminate called after throwing an instance of 'std::bad_alloc'
what():  std::bad_alloc
Aborted (core dumped)

```

May there can be an issue in reference of wheat that I am using, I got this reference from here. Here are the first ten lines of the reference using head command:

```
> head wheat.fa
>Pt dna:chromosome chromosome:TGACv1:Pt:1:135835:1 REF
ACAGAAATACCCAATATCTTGTTCTAGCAAGATATTGGGTATTTCTGTCTTTTCTTTCTT
CAAAAATTCTTATATGTTAGCAGAAAAACCTTATCCATTAATAGATGGAACTTCAACAGC
AGCTAAGTCTAGAGGGAAGTTGTGAGCATTACGTTCGTGCATTACTTCCATACCAAGGTT
AGCACGGTTGATGATATCAGCCCAAGTATTAATAACGCGACCTTGACTATCAACTACAGA
TTGGTTGAAATTGAAACCATTTAGGTTGAAAGCCATAGTACTAATACCTAAAGCAGTGAA
CCAGATTCCTACTACAGGCCAAGCAGCCAAGAAGAAGTGTAAAGAACGAGAGTTGTTGAA
ACTAGCATATTGGAAGATTAATCGGCCAAAATAACCATGAGCAGCCACAATATTATAAGT
TTCTTCCTCTTGACCAAATTTGTAACCCTCATTAGCAGATTCATTTTCAGTAGTTTCCCT
GATCAAACTAGAGGTTACCAAGGAACCATGCATAGCACTGAATAGGGAACCGCCGAATAC

```

How this problem can be resolved?

# Answer
> 3 votes
In addition to the other answer (to reduce the number of threads), you might also need to add a `--genomeChrBinNbits` argument. Your reference has more than 5000 scaffolds, so you need to adjust the `--genomeChrBinNbits` argument to reduce RAM consumption. Found here.

They recommend:

> --genomeChrBinNbits} = min(18, log2(GenomeLength/NumberOfReferences))

In your case:

Genome size 17G, number of scaffolds 735945

So `--genomeChrBinNbits` should be set to log2(17,000,000,000/735,945) = 14.4955, so let's say 14 or 15.

# Answer
> 5 votes
Try reducing the number of threads.

When multithreading the index creation (and other memory bound tasks), the memory usage increases linearly with the number of threads. If you required a 10 GiB working space work with a single thread, then the requirement for 10 threads would be at 100 GiB.

STAR is particularly memory hungry for the index creation to start with, and your genome is exceptionally large. I’ve had STAR require \> 32 GiB of memory for the human genome, with a moderate number of threads. I’m not surprised that it fails with a 4× larger genome and many more threads: it would probably require at least 2× as much RAM as your machine has.

---
Tags: rna-seq, rna-alignment, star
---


title: Preparing data for UpSet Plot in R

text:
# Question
Title: Preparing data for UpSet Plot in R
As i have multiple comparison so I m trying this `UpSet` plot, but it needs data to be formatted in particular way, one of the ways is as a list, so now the issue is since between one comparison the no of genes coming out may not be same as the other one, so there is a difference in dimension so if i bind them as data frame then i get dimension error, so the best way to get through is make list so what i have tried so far

```
A <- read.csv("GMP_MONO_UP_DOWN_CDS_EPI.txt",header = TRUE)
B <- read.csv("CMP_GMP_UP_DOWN_EPI.txt",header = TRUE)
C <- read.csv("HSC_BLAST_UPDOWN_EPI.txt",header = TRUE,sep = '\t')

dim(A)
dim(B)
ONE <- A[c(1)]#############just taking out the gene column########
TWO <- B[c(1)]#############just taking out the gene column########
THREE <- B[c(1)]#############just taking out the gene column########
dim(ONE)
dim(TWO)
b1 <- apply(ONE,2,list)
b2 <- apply(TWO,2,list)

```

Now this is the code from their help page

```
upset(movies, nsets = 7, nintersects = 30, mb.ratio = c(0.5, 0.5),
      order.by = c("freq", "degree"), decreasing = c(TRUE,FALSE))

```

so the data "movies" class is data frame, there is not dimension issue as such as in my case if i want to intersect two comparison which may not have the same length, so how do i overcome the issue any suggestion or help would be highly appreciated. Update my question with one more data to add.

# Answer
> 3 votes
Don't you just need a binary matrix with presence/absence data for all of the genes?

```
library(dplyr)
combined <- full_join(data.frame(gene=ONE, GMP_MONO=1),
                      data.frame(gene=TWO, CMP=1))
combined[is.na(combined)] <- 0

```

Edit: To loop over a large number of data frames, the easiest solution is to use the `reduce` function from the purrr package:

```
combined <- reduce(list(data.frame(gene=c('gene1', 'gene2', 'gene3'), set1=1),
                        data.frame(gene=c('gene3', 'gene4', 'gene5'), set2=1),
                        data.frame(gene=c('gene1', 'gene4'), set3=1)
                       ), full_join)
combined[is.na(combined)] <- 0
combined
#    gene set1 set2 set3
# 1 gene1    1    0    1
# 2 gene2    1    0    0
# 3 gene3    1    1    0
# 4 gene4    0    1    1
# 5 gene5    0    1    0

```

# Answer
> 7 votes
We can plot using lists as well, using `fromList` function, see example:

```
library(UpSetR)

myGeneSets <- list(
  set1 = c("gene1","gene2","gene3"),
  set2 = c("gene1","gene4","gene5", "gene7"),
  set3 = c("gene1","gene7")
  )

upset(fromList(myGeneSets))

```

> **fromList**
>
> A function to convert a list of named vectors to a data frame compatible with UpSetR.

```
fromList(myGeneSets)
#   set1 set2 set3
# 1    1    1    1
# 2    1    0    0
# 3    1    0    0
# 4    0    1    0
# 5    0    1    0
# 6    0    1    1

```

# Answer
> 1 votes
You need to set a common number of rows, one for each gene (make sure that the gene identifiers are the same for both files). Then in the first column you can set 1 for the file A, and 0 otherwise. In the second column you can do the same. Then you add this dataframe as the input of the `upset` function.

```
genes <- unique(c(ONE, TWO))
df <- data.frame(A = ifelse(genes %in% A, 1, 0), B = ifelse(genes %in% b, 1, 0))
upset(df, ...)

```

---
Tags: r, upsetr
---


# Question
Title: How to subset a GRanges via an argument passed into a function?
Let's say I have following example GRanges:

```
> library(GenomicRanges)
> gr = GRanges(
+     seqnames = Rle(c("chr1", "chr2", "chr1", "chr3"), c(1, 3, 2, 4)),
+     ranges = IRanges(101:110, end = 111:120, names = head(letters, 10)),
+     strand = Rle(strand(c("-", "+", "*", "+", "-")), c(1, 2, 2, 3, 2)),
+     score = 1:10,
+     GC = seq(1, 0, length=10))
> gr
GRanges object with 10 ranges and 2 metadata columns:
    seqnames     ranges strand |     score                GC
       <Rle>  <IRanges>  <Rle> | <integer>         <numeric>
  a     chr1 [101, 111]      - |         1                 1
  b     chr2 [102, 112]      + |         2 0.888888888888889
  c     chr2 [103, 113]      + |         3 0.777777777777778
  d     chr2 [104, 114]      * |         4 0.666666666666667
  e     chr1 [105, 115]      * |         5 0.555555555555556
  f     chr1 [106, 116]      + |         6 0.444444444444444
  g     chr3 [107, 117]      + |         7 0.333333333333333
  h     chr3 [108, 118]      + |         8 0.222222222222222
  i     chr3 [109, 119]      - |         9 0.111111111111111
  j     chr3 [110, 120]      - |        10                 0
  -------
  seqinfo: 3 sequences from an unspecified genome; no seqlengths

```

It is straightforward to subset this GRanges as one would subset an R data.frame, e.g.

```
> gr$score
 [1]  1  2  3  4  5  6  7  8  9 10

```

However, I am writing a function where the user would pass through a string for the metadata column to subset e.g.

```
subsetter = function(granges, column){
    return(granges$column)
}

```

I try to use the function as follows:

```
> subsetter(gr, score)

```

which outputs `NULL`

This obvious doesn't work, as it will subset the `gr` by `column`, not `score`.

```
> subsetter(gr, score)
NULL

```

If this was a data.frame, I could subset as follows:

```
gr[[(score)]]

```

but this gives the error

```
Error in gr[[(score)]] : this S4 class is not subsettable

```

Are there any other ways to subset a GRanges besides using `$`?

# Answer
> 6 votes
```
subsetter = function(gr, cname) {
    return(mcols(gr)[[cname]])
}

```

You can then use things like `subsetter(gr, "GC")` and `subsetter(gr, "score")`.

# Answer
> 5 votes
In such cases it is best to check the `str` of an object *gr*. Then we could see that meta data is just a *dataframe* inside *S4* object:

```
gr@elementMetadata
# DataFrame with 10 rows and 2 columns
#        score        GC
#    <integer> <numeric>
# 1          1 1.0000000
# 2          2 0.8888889
# ...

```

Once we know how to access the *dataframe*, we can then subset the same way as we would subset a *dataframe*:

```
gr@elementMetadata[, "score" ]

```

I don't see the point of wrapping this into a function, but if it is needed:

```
subsetter <- function(gr, cname) gr@elementMetadata[, cname ]

# test the function
subsetter(gr, "score")
# [1]  1  2  3  4  5  6  7  8  9 10

```

---

**Edit**

Alternatively, we could use `GenomicRanges::elementMetadata()` function to extract *dataframe* then subset:

```
subsetter <- function(gr, cname) elementMetadata(gr)[, cname ]

```

# Answer
> 4 votes
In addition to Devon’s answer, you can use non-standard evaluation to make this work with R symbol names rather than string arguments:

```
subsetter = function (gr, cname) {
    cname = as.character(substitute(cname))
    mcols(gr)[[cname]]
}

```

Now using `subsetter(gr, score)` works analogously to `gr$score`, which may be desirable. But, likewise, you’ll no longer be able to use a variable:

```
var = 'score'
subsetter(gr, var) # fails, same as `gr$var`

```

---
Tags: r, bioconductor, genomicranges
---