title: Tree cut issue in WGCNA

text:
# Question
Title: Tree cut issue in WGCNA
```
library(WGCNA);
# The following setting is important, do not omit.
library(flashClust)
enableWGCNAThreads()
options(stringsAsFactors = FALSE);
#Read in the female liver data set
femData = read.csv("CTRL_WGCNA.txt",sep = "\t");

head(femData)
# Read in the male liver data set
maleData = read.csv("BPD_WGCNA.txt",sep = "\t");
# Take a quick look at what is in the data sets (caution, longish output):
dim(femData)
names(femData)
dim(maleData)
names(maleData)

#=====================================================================================
#
#  Code chunk 2
#
#=====================================================================================

# We work with two sets:
nSets = 2;
# For easier labeling of plots, create a vector holding descriptive names of the two sets.
setLabels = c("Female liver", "Male liver")
shortLabels = c("Female", "Male")
# Form multi-set expression data: columns starting from 9 contain actual expression data.
multiExpr = vector(mode = "list", length = nSets)

multiExpr[[1]] = list(data = as.data.frame(t(femData[-c(1)])));
names(multiExpr[[1]]$data) = femData$Gene;
rownames(multiExpr[[1]]$data) = names(femData)[-c(1)];
multiExpr[[2]] = list(data = as.data.frame(t(maleData[-c(1)])));
names(multiExpr[[2]]$data) = maleData$Gene;
rownames(multiExpr[[2]]$data) = names(maleData)[-c(1)];
# Check that the data has the correct format for many functions operating on multiple sets:
exprSize = checkSets(multiExpr)

#=====================================================================================
#
#  Code chunk 3
#
#=====================================================================================

# Check that all genes and samples have sufficiently low numbers of missing values.
gsg = goodSamplesGenesMS(multiExpr, verbose = 3);
gsg$allOK

#=====================================================================================
#
#  Code chunk 4
#
#=====================================================================================

if (!gsg$allOK)
{
  # Print information about the removed genes:
  if (sum(!gsg$goodGenes) > 0)
    printFlush(paste("Removing genes:", paste(names(multiExpr[[1]]$data)[!gsg$goodGenes],
                                              collapse = ", ")))
  for (set in 1:exprSize$nSets)
  {
    if (sum(!gsg$goodSamples[[set]]))
      printFlush(paste("In set", setLabels[set], "removing samples",
                       paste(rownames(multiExpr[[set]]$data)[!gsg$goodSamples[[set]]], collapse = ", ")))
    # Remove the offending genes and samples
    multiExpr[[set]]$data = multiExpr[[set]]$data[gsg$goodSamples[[set]], gsg$goodGenes];
  }
  # Update exprSize
  exprSize = checkSets(multiExpr)
}

#=====================================================================================
#
#  Code chunk 5
#
#=====================================================================================

sampleTrees = list()
for (set in 1:nSets)
{
  sampleTrees[[set]] = hclust(dist(multiExpr[[set]]$data), method = "average")
}

#=====================================================================================
#
#  Code chunk 6
#
#=====================================================================================

pdf(file = "SampleClustering.pdf", width = 12, height = 12);
par(mfrow=c(2,1))
par(mar = c(0, 4, 2, 0))
for (set in 1:nSets)
  plot(sampleTrees[[set]], main = paste("Sample clustering on all genes in", setLabels[set]),
       xlab="", sub="", cex = 0.7);
dev.off();

  #=====================================================================================
#
#  Code chunk 7
#
#=====================================================================================

# Choose the "base" cut height for the female data set
baseHeight = 80
# Adjust the cut height for the male data set for the number of samples
cutHeights = c(80, 100*exprSize$nSamples[2]/exprSize$nSamples[1]);
# Re-plot the dendrograms including the cut lines
pdf(file = "SampleClustering.pdf", width = 12, height = 12);
par(mfrow=c(2,1))
par(mar = c(0, 4, 2, 0))
for (set in 1:nSets)
{
  plot(sampleTrees[[set]], main = paste("Sample clustering on all genes in", setLabels[set]),
       xlab="", sub="", cex = 0.7);
  abline(h=cutHeights[set], col = "red");
}
dev.off();

  #=====================================================================================
#
#  Code chunk 8
#
#=====================================================================================

for (set in 1:nSets)
{
  # Find clusters cut by the line
  labels = cutreeStatic(sampleTrees[[set]], cutHeight = cutHeights[set])
  # Keep the largest one (labeled by the number 1)
  keep = (labels==1)
  multiExpr[[set]]$data = multiExpr[[set]]$data[keep, ]
}
collectGarbage();
 # Check the size of the leftover data
exprSize = checkSets(multiExpr)
exprSize

$nSets
[1] 2

$nGenes
[1] 8213

$nSamples
[1] 0 0

$structureOK
[1] TRUE

```

I think the issue is with the tree cut height because when i check the label all of them are zeroes i guess that is issue so far which is why im getting exprSize nsamples as `0 0`instead of my actual sample size.

I have tried it with different cut height but so far no success may its a pretty simple issue but im not able to figure it out.

Data which im using here is the link

Any suggestion or help would be highly appreciated

# Answer
It is likely that the cut heights set in the code (80 and 100\*ratio of sample numbers) are too low. Look at the sample trees plotted into a pdf file (SampleClustering.pdf) and note the merging heights. You need to set the cut heights appropriately so that most samples go into a single cluster.

Peter
> 2 votes

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


title: read and matching pattern with python

text:
# Question
Title: read and matching pattern with python
I have read the content of a text file into pandas and needed some help matching the pattern.

Here is the pattern, where `a` can any number greater than zero and `b` can be zero or any number greater than zero:

```
father=a|0 || 0|a
mother=a|0 || 0|a
daughter1=b|0 || 0|b
daughter2=a|0 || 0|a
daughter3=b|0 || 0|b
son1=a|0 || 0|a
son2=b|0 || 0|b

```

Here is the function that reads the content of text file to pandas

```
import pandas as pd
def read():
    varian=pd.read_csv('file.txt',sep = "\t")

father=varian["FATHER"]
mother=varian["MOTHER"]
daughter1=varian["DAUGHTER1"]
daughter2=varian["DAUGHTER2"]
daughter3=varian["DAUGHTER3"]
son1=varian["SON1"]
son2=varian["SON2"]

```

The content of the file looks like this

```
HROM    POS REF ALT FATHER  MOTHER  DAUGHTER1   DAUGHTER2   DAUGHTER3   SON1    SON2    INFO    FREQUECY
1   1226852 G   C   1/0 0/0 1|0 1|0 1|0 1|0 0|0 AN=2184;AC=12   0.005494505
1   1847936 C   T   0/1 1/1 1|1 1|1 1|1 1|1 0|1 AC=1;AN=2184    0.000457875
1   2428427 C   G   0/1 0/0 0|0 0|0 0|0 0|0 1|0 AC=4;AN=2184    0.001831502
1   2515616 G   A   1/0 1/0 1/0 1|1 1|1 1/0 1/0 AC=1;AN=2184    0.000457875
1   3801895 T   C   1/0 1/0 0|0 1/0 1/0 0|0 1|1 AC=10;AN=2184   0.004578755

```

# Answer
> 1 votes
This is not a bioinformatics question, but a basic pandas question, so SO is a better place for this... But your odd looking VCF file can be parsed as follows without even using regex:

```
import pandas as pd

def splitter(cell, key):
    parts = cell.split('||')
    if key in parts[0]:
         p = parts[0].split('|')
    else:
         p = parts[1].split('|')
    sp = [x.strip() for x in p if key not in x]
    return float(sp[0])

def read():
    splitter_a = lambda v: v.split('||')[0].split('|')[1])
    varian=pd.read_csv('file.txt',sep = "\t")
    a_splitter = lambda c: splitter(c, 'a')
    b_splitter = lambda c: splitter(c, 'b')
    varian = varian.assign(father_a=varian.FATHER.apply(a_splitter)\
              .assign(father_a=varian.FATHER.apply(b_splitter)\
              .assign(mother_a=varian.MOTHER.apply(la_splitter)\

```

etc.

`display(varian)` or `print(varian.columns)` will show that you have new columns called `father_a` etc.

# Answer
> 1 votes
Its easy. You want to remove homozygotes and here is the solution for the FATHER class. Note the MOTHER class will comprise homozygotes in the answer because you are only searching for FATHER. To perform an intersection one approach is just to concatenate the mask.

```
import pandas as pd

genedict = {\
            'FATHER': ['1/0', '1/1', '0/0', '1/0', '0/1'], \
            'MOTHER': ['1/1', '1/1', '1/0', '1/0', '0/0']}

df = pd.DataFrame(genedict)
mask = df.FATHER.str.contains(r'1[|/]0') + df.FATHER.str.contains(r'0[|/]1')
df[mask]

```

**OUTPUT**

```
  FATHER MOTHER
0    1/0    1/1
3    1/0    1/0
4    0/1    0/0

```

Note, this is not your original dataset, its a dictionary to demonstrate the principle and note all your other columns will carry across unless you want to 'drop' them. Finally, pandas has its own regex system, regardless the traditional r'' code is still needed

---
Tags: python
---


# Question
Title: What value of reference RMSD after docking simulation from any starting coordinate would be considered outrageous to have?
I need to lightly verify a set of docking simulations from the reference RMSD of the outputs, and as we all know reference RMSD depends on the initial coordinates of the ligand set. What value of reference RMSD after docking from any starting coordinate would be considered outrageous to have?

# Answer
Let's get some definitions out of the way:

**RMSD**, root mean square deviation, is a metric of distance of molecule A and molecule B. Think of it as average euclidean distance. It us casually called just 'deviation', but has nothing to do with standard deviation (the root of variance, 2nd moment). The units are Ångström.

**Reference RMSD** is a term that appears in some manuals and is a the RMSD between the calculated pose and the starting pose. When you set up your docking you as a user define the starting position and how much it can deviate from that position. If you have a crystal structure with a ligand and you are docking the same or similar compound, then yes, it makes total sense to worry about this metric. **In which case, RMSD greater than 2 are worrisome in a 2 Å resolution structure.** You especially see RMSD as a metric of accuracy when PDB structures are redocked, wherein the same ligand is docked and reference is the file from the PDB (or its PDB Redo derivative), so RMSD greater than 2 Å is bad as it would not match the density, while under 1 is really good —note that in tests about 50% and 75% of re-docked ligands are docked less than 2 with AutoDock (free) and ICM-Dock (pay) respectively. But if you are doing it for an experiment and not for a test, you would probably have set the docking to be highly constrained in the first place (the grid size in AutoDock or position randomisation factor in most others, eg. Rosetta Ligand dock). If you actually care about reference RMSD the thing to manually look at is the top 3-4 poses in the electron density of your native structure. If you don't really know how your ligand docks please ignore this value as it is just a product of your starting pose. It is *not* a measure of how "happy" the complex is.

**ΔGibbs** is a negative value in kcal/mol that is basically a metric of how "happy" the atoms are. It is actually an energy potential, but confer Wikipedia pages for the proper definition. The more negative it is the happier the atoms are. Force-field calculations solve this value and the difference between undocked and docked (ΔΔG) is a good metric of how well your ligand docks. Please use this.

**Affinity** ΔΔG is not a k\_off or a k\_on. That is a result of the ligand binding to nearby positions and how well it can get into the binding site (*e.g.* in P450 enzymes the ligand orientation is dictated by how it squeezes into the entrance way). And how insoluble the ligand is, which is a major problem\*. For a better metric, there is **dynamic undocking** where a constrained MD simulation is done with the ligand getting pulled out of the pocket. But this is not trivial.

∗ As any comp-chemist can tell you it is mighty hard to calculate solubility. Some force-field calculations do not even take account of π-π or S-π interactions. Drug-like compounds are rigid aromatic and rigid and highly hydrophobic. Many programs use a hack called implicit waters, where the water is not simulated as individual molecules but as a flowing thing. This means that ΔΔG benefits of constrained waters are ignored and alternative protontation states of water too. This is fine for protein (mostly), but has major errors for drugs. But in terms of thing to look out for in a docking algorithm, explicit waters is not the first as conformer sampling and side chain perturbation are more important.
> 2 votes

# Answer
2A would be huge for docking, but agreed depends on the resolution.

Its 10A if you are looking at protein-protein bonds in a beta pleated sheet and your molecules need to be a lot smaller than that.
> 1 votes

---
Tags: proteins, docking
---


title: WGCNA module preservation analysis

text:
# Question
Title: WGCNA module preservation analysis
I have case control studies,4 different condition where I have 1 control and against which I'm have 3 comparisons to make.

I did WGCNA for individuals samples .Here I have the kept the number of genes equal in all the cases .

Now my question is whether the number of modules which are obtained from WGCNA would be similar in number if I run module preservation between two condition .

Lets say i get 15 modules in Control condition and 16 treated condition if I run module preservation the number of modules which would come would be 15+16= 31 or its independent .

The power I have kept same for both WGCNA and module preservation analysis. But Still I get different module numbers.

Any suggestion would be highly appreciated

# Answer
> 2 votes
Your question is not clear, you seem to mix the relevant terminology...

First, you cannot run WGCNA for individual samples, only for individual conditions (with or without WT, depending on what question you're trying to answer and assuming you have enough samples in each condition). If that's what you have done, fine. The number of modules in each independent analysis will generally be different, that's because each condition is different. The module identification built into WGCNA as the default approach does not allow the user to specify the number of modules beforehand, but if you insist on getting a pre-specified number of modules, you could run WGCNA using the "manual" way (for example, Tutorial I, section 2b) and substitute another clustering method (e.g., PAM) that allows you to specify the number of clusters. You may find the results inferior though.

Second, module preservation does not result in new modules. Module preservation evaluates whether modules found in data set A are preserved in a different data set B. To get new, common modules, you could run a consensus network analysis on pairs of conditions or all 3 conditions together. The number of modules you would find in a consensus analysis cannot be predicted from the number of modules you find in individual data sets.

---
Tags: wgcna
---


# Question
Title: What do the symbols mean in minimap2's gap cost equation?
Heng Li's paper on minimap2 is understandable to me right up until the moment that it gives this definition of gap cost:

> ### 2.2.1 Alignment with 2-piece affine gap cost
>
> Minimap2 performs DP-based global alignment between adjacent anchors in a chain. It uses a 2-piece affine gap cost (Gotoh, 1990):
>
> $$\gamma\_a(l) = \min \\{ q + |l| \cdot e, \tilde{q} + |l| \cdot \tilde{e} \\}$$

Suddenly, out of nowhere, we've got *five* variables used without prior introduction: $l$, $q$, $e$, $\tilde{q}$, and $\tilde{e}$.

Heng doesn't state their meanings, and based on a quick skim and text search of the cited paper by Gotoh, that paper doesn't use these variable names anywhere either.

What does this equation mean?

# Answer
> 5 votes
Here is the translation:

* $l$ = gap length
* $q$ = initial gap penalty for *shorter* gaps
* $e$ = *short* gap extension penalty per every additional gap
* $\tilde{q}$ = initial gap penalty for *longer* gaps
* $\tilde{e}$ = *long* gap extension penalty per every additional gap

To understand the first of the two expressions being minimized, see the Wikipedia section on affine gap penalties.

Notice the similar notation structure

$$A + B \cdot L$$

where

* $A$ is the known gap penalty,
* $B$ is the gap extension penalty, and
* $L$ is the length of the gap

To understand the other variables with the tildes on top, see the text below Equation 2 (emphasis added to key explanation):

> On the condition that $e \> \tilde{e}$, it applies cost $q + |l| \cdot e$ to gaps shorter than $\[(\tilde{q} - q)/(e - \tilde{e})\]$ and applies $\tilde{q} + |l| \cdot \tilde{e}$ to longer gaps. **This scheme helps to recover longer insertions and deletions (INDELs)**.

In other words, we make the assertion that we should penalize longer gaps differently than shorter gaps.

As a bonus, here's a visual understanding of why we use the minimization function on those two expressions using Python and some made-up (and probably not biologically useful) penalty terms. **Note** the requirement that $e \> \tilde{e}$.

```
# Import libraries
import numpy as np
import matplotlib.pyplot as plt

# Setup
l = np.array(range(0, 15, 1))
y_short = eval('10 + 5*l') # Q=10, E=5
y_long = eval('20 + 3*l')  # Q=20, E=3

# Plot
plt.plot(l, y_short)
plt.plot(l, y_long)
plt.legend(["Short Gap Penalty", "Long Gap Penalty"])

```

Using our example and the expression $\[(\tilde{q} - q)/(e - \tilde{e})\]$, it will apply the first expression in the minimization function (the one without the tildes) when the gap length is shorter than 5.

$$ \[(\tilde{q} - q)/(e - \tilde{e})\] = \[20 - 10)/(5 - 3)\] = 10 / 2 = 5 $$

We can see from our plot that the two lines intersect at 5, which we've verified through some arithmetic.

---
Tags: sequence-alignment, read-mapping, minimap2
---


title: How I deal with this expression set?

text:
# Question
Title: How I deal with this expression set?
I have a gene expression raw counts like below for 16 patients

```
         A  B    C      D   E   F
A2M     511 1623 665    208 553 469
AADAT   34  137 372     7   52  124
ABCB1   119 114 123    22   25  89
ABCB11  27  186 200    27   30  49
ABCC2   1   10  21     1    8   3

```

For each patients I also have some characteristics

```
Patients                          IC          TCs
A1                                3            3
B1                                2            0
C1                                2            0
D1                                0            0

```

I want to find genes differential expressed in these patients regarding these characteristics but I don't know how to do that

Any help please?

# Answer
> 1 votes
Look up a tutorial on DESeq2. To keep things simple for your first time, I'd reduce your phenotypes to "high" and "low" bins instead of trying to use the cell count numbers as they are. Then once you understand what the software is doing, you can get fancier.

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


title: fixing chromosome with new contig

text:
# Question
Title: fixing chromosome with new contig
I want to fix a chromosome of around 50-100Mb (`chr.fasta`) for a region that is missing a few kilobases of data. I have the missing sequence with a window upstream and downstream that matches perfectly to the reference (`contig.fasta`).

What's the best way to produce a new chromosome, `fixed.fasta` that contains the new data from `contig.fasta`?

See depiction below:

```
chr.fasta    |---------------------------------------------------   -------------|
contig.fasta                                                   |-===-|

```

Final intended result:

```
fixed.fasta  |---------------------------------------------------===-------------|

```

Any ideas?

# Answer
> 1 votes
I have not received any answer to the question, so I created a small tool to fix the file myself.

I am assuming as per the question that there are perfect matches upstream and downstream to the new sequence, and that anything in between can be deleted/updated by the new sequence.

Given the first\_n and last\_n bases of the diagram in the question, and the entire patch sequence, we can create a single-line string of the reference with `seqtk seq -l9999999 -` and then edit it with the line below (in Perl regex):

```
$chromosome_str =~ s/$first_n.+$last_n/$patch/i

```

This creates the new string that can be used as the new fasta file.

# Answer
> 0 votes
You can give a try with `megamerger` from the `emboss` suite if you have (almost) perfect sequence match.

In case the gap within your `chr.fasta` complicates the sequence alignment and hence the resulting contig, you can generate a "tiling path" with: upstream of the gap from `chr.fasta` \+ `contig.fasta` \+ downstream of the gap from `chr.fasta` (there should be sequence overlaps of course)

In a (more realistic) scenario where there is no perfect sequence match between the overlapping regions, you can fall back to `emboss merger`, however, it is not designed to work on large sequences.

---
Tags: assembly, chromosomes
---


# Question
Title: Is there an easy way to create a summary of a VCF file (v4.1) with structural variations?
I got a bunch of vcf files (v4.1) with structural variations of bunch of non-model organisms (i.e. there are no known variants). I found there are quite a some tools to manipulate vcf files like VCFtools, R package vcfR or python library PyVCF. However none of them seems to provide a quick summary, something like (preferably categorised by size as well):

```
type    count
DEL     x
INS     y
INV     z
....

```

Is there any tool or a function I overlooked that produces summaries of this style?

I know that vcf file is just a plain text file and if I will dissect `REF` and `ALT` columns I should be able to write a script that will do the job, but I hoped that I could avoid to write my own parser.

--- edit ---

So far it seems that only tool that aims to do summaries (@gringer answer) is not working on vcf v4.1. Other tools would provide just partial solution by filtering certain variant type. Therefore I accept my own parser perl/R solutions, till there will be a **working** tool for **stats** of vcf with **structural variants**.

# Answer
The "my own parser" solutions. The information I was searching for in part of column `INFO`, namely in variables `SVLEN` and `SVTYPE`.

very fast SV types + counts (by @user172818 in comment) :

```
zcat var.vcf.gz | perl -ne 'print "$1\n" if /[;\t]SVTYPE=([^;\t]+)/' | sort | uniq -c

```

quite slow SV types + counts + sizes :

```
SV_colnames <- c('CHROM', 'POS', 'ID', 'REF', 'ALT', 'QUAL', 'FILTER', 'INFO', 'FORMAT', 'SAMPLE1')

ssplit <- function(s, split = '='){
    unlist(strsplit(s, split = split))
}

# note, capital letters just respect original naming conventions of the VCF file
getSVTYPE <- function(info){
    ssplit(grep("SVTYPE", info, value = T))[2]
}

getSVLEN <- function(info){
    SVLEN <- ssplit(grep("SVLEN", info, value = T))
    ifelse(length(SVLEN) == 0, NA, as.numeric(SVLEN[2]))
}

load_sv <- function(file){
    vcf_sv_table <- read.table(file, stringsAsFactors = F)
    colnames(vcf_sv_table) <- SV_colnames
    # possible filtering
    # vcf_sv_table <- vcf_sv_table[vcf_sv_table$FILTER == 'PASS',]
    vcf_sv_table_info <- strsplit(vcf_sv_table$INFO, ';')
    vcf_sv_table$SVTYPE <- unlist(lapply(vcf_sv_table_info, getSVTYPE))
    vcf_sv_table$SVLEN <- unlist(lapply(vcf_sv_table_info, getSVLEN))
    return(vcf_sv_table)
}

sv_df <- load_sv('my_sv_calls.vcf')
table(sv_df$SVTYPE)

```
> 6 votes

# Answer
According to the `bcftools` man page, it is able to produce statistics using the command `bcftools stats`. Running this myself, the statistics look like what you're asking for:

```
# This file was produced by bcftools stats (1.2-187-g1a55e45+htslib-1.2.1-256-ga356746) and can be plotted using plot-vcfstats.
# The command line was: bcftools stats  OVLNormalised_STARout_KCCG_called.vcf.gz
#
# Definition of sets:
# ID    [2]id   [3]tab-separated file names
ID  0   OVLNormalised_STARout_KCCG_called.vcf.gz
# SN, Summary numbers:
# SN    [2]id   [3]key  [4]value
SN  0   number of samples:  108
SN  0   number of records:  333
SN  0   number of no-ALTs:  0
SN  0   number of SNPs: 313
SN  0   number of MNPs: 0
SN  0   number of indels:   20
SN  0   number of others:   0
SN  0   number of multiallelic sites:   0
SN  0   number of multiallelic SNP sites:   0
# TSTV, transitions/transversions:
# TSTV  [2]id   [3]ts   [4]tv   [5]ts/tv    [6]ts (1st ALT) [7]tv (1st ALT) [8]ts/tv (1st ALT)
TSTV    0   302 11  27.45   302 11  27.45
# SiS, Singleton stats:
...
# IDD, InDel distribution:
# IDD   [2]id   [3]length (deletions negative)  [4]count
IDD 0   -9  1
IDD 0   -2  4
IDD 0   -1  6
IDD 0   1   4
IDD 0   2   1
IDD 0   3   3
IDD 0   4   1
# ST, Substitution types:
# ST    [2]id   [3]type [4]count
ST  0   A>C 2
ST  0   A>G 78
ST  0   A>T 2
ST  0   C>A 5
ST  0   C>G 0
ST  0   C>T 66
ST  0   G>A 67
ST  0   G>C 0
ST  0   G>T 1
ST  0   T>A 1
ST  0   T>C 91
ST  0   T>G 0
# DP, Depth distribution
# DP    [2]id   [3]bin  [4]number of genotypes  [5]fraction of genotypes (%)    [6]number of sites  [7]fraction of sites (%)
DP  0   >500    0   0.000000    333 100.000000

```

This is annotating what is explicitly in the VCF file, and that's SNVs (and INDELs). If you want a structural variant analysis (i.e. on a larger scale than single nucleotides), then you'll need to use something that does more than a summary of the VCF file.

While Inversions, large-scale deletions, and breakpoints are part of the 4.1 specifications, they are unfortunately not currently supported by bcftools.
> 9 votes

# Answer
You might want to try Bio-VCF. From the authors description

> Bio-vcf is a new generation VCF parser, filter and converter. Bio-vcf is not only very fast for genome-wide (WGS) data, it also comes with a really nice filtering, evaluation and rewrite language and it can output any type of textual data, including VCF header and contents in RDF and JSON.

In addition it is a very fast parser. The rewrite DSL might help you customise your filtering and needs.
> 2 votes

# Answer
For structural variants I think you can use SURVIVOR like `SURVIVOR stats` which was design specifically for such purpose (statistics on SV file).
> 2 votes

# Answer
Let's say you have an INFO field called `SVTYPE` to indicate the structure variant type.

Here is how you get the stats easily:

1. Install `vcfstats`: `pip install vcfstats`
2. Define a macro to extract the `SVTYPE` information, say in `svtype.py`:

```
from vcfstats.macros import cat
@cat
def SVTYPE(variant):
    return variant.INFO['SVTYPE']

```
3. Generate the statistics:

```
vcfstats --vcf <your vcf file> \
   --outdir <the output directory> \
   --macro svtype.py \
   --formula 'COUNT(1) ~ SVTYPE[DEL,INS,INV...]' \ # desired svtypes
   --title 'Number of SV types'

```

Then in the output directory you will have the stats file named `Number_of_SV_Types.txt` and a plot for it as well.

Check the details at https://github.com/pwwang/vcfstats
> 1 votes

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


title: BED file from .bam alignment structure

text:
# Question
Title: BED file from .bam alignment structure
I am using a pipeline from nanopore to detect structural variants (SVs) in a human sample with long-reads sequencing. The first steps of the pipeline are:

1. index the reference genome with `minimap2`
2. map the long reads to a reference genome using same tool as above
3. Get BED coordinates from the `.bam` file.

While the first two steps uses `minimap`, the third one uses a python script written by the developers of the pipeline, which can be found here.

After the pipeline calls `bamref2bed` i.e. the python script above, my `target.bed` looks like the following:

```
chr1    0   249250621
chr2    0   243199373
chr3    0   198022430
chr4    0   191154276
chr5    0   180915260
chr6    0   171115067
chr7    0   159138663
chr8    0   146364022
chr9    0   141213431
chr10   0   135534747
chr11   0   135006516
chr12   0   133851895
chr13   0   115169878
chr14   0   107349540
chr15   0   102531392
chr16   0   90354753
chr17   0   81195210
chr18   0   78077248
chr19   0   59128983
chr20   0   63025520
chr21   0   48129895
chr22   0   51304566
chrX    0   155270560
chrY    0   59373566

```

As far as I am concerned the first three columns of a bed file refer to (I) chr name, (II) start position (III) end position. The main thing I am wondering is how come all chromosomes have a start position of 0. Also, reading about converting bam to bed files I guess the most important feature would be the coverage across the whole chromosomes, which in this file is missing.

But my main question would be: is it normal to have a start position of 0 for all chromosomes or is there something wrong here?

# Answer
> 4 votes
I cannot access the pipeline from nanopore, but looking at the python script you provided:

```
with pysam.AlignmentFile(bam_file, "r") as bam:
    for tid in range(0, bam.nreferences):
        ref_name = bam.get_reference_name(tid)
        if len(list(filter(lambda x: x in ref_name, filter_names))) > 0:
            continue

        print(ref_name,
              "0",
              str(bam.get_reference_length(ref_name)),
              sep='\t')

```

So the script gets the length of the chromosome and by default, outputs the start as 0.Zero based coordinate system is sometimes used in UCSC systems, so I am only guessing the authors of the pipeline use this convention somewhere further down.

As for the coverage you mentioned, that is usually output as a score in the 5th column (bed format). In those situations, usually the coordinates of the alignments in the bam file are converted into bed format. In this pipeline, i think the authors are simply extracting the chromosome information out and putting it into a tabular form.

---
Tags: sequence-alignment, nanopore, bed, long-reads, structural-variation
---


title: RNAseq data analysis

text:
# Question
Title: RNAseq data analysis
I am currently doing a Differential Gene Expression of RNAseq data in Lung Adenocarcinoma using TCGAbiolinks. In the data Preprocessing step TCGAanalyze\_Normalisation, I am confused as to which method to select for normalization (geneLength or gcContent). What are the possible implications of both methods and how to decide which is the best suited model for the analysis?

Is there a better way to perform the DGE ?

# Answer
The TCGAbiolinks reference manual provides a clear explanation of the two methods:

> Normalization for RNA-Seq Numerical and graphical summaries of RNA-Seq read data. Withinlane normalization procedures to adjust for GC-content effect (or other gene-level effects) on read counts: loess robust local regression, global-scaling, and full-quantile normalization (Risso et al., 2011). Between-lane normalization procedures to adjust for distributional differences between lanes (e.g., sequencing depth): global-scaling and full-quantile normalization (Bullard et al., 2010).

Ideally, you would know if your samples spanned different lanes or not, which would indicate which you should use - "geneLength" if so, "gcContent" if not. If you're using TCGA data, that metadata is likely available somewhere.
> 2 votes

---
Tags: rna-seq
---


# Question
Title: Tag density plot form chip or atac seq data
I m quoting this paper , figure B where they write

"B) Violin plots for the normalized ATAC-seq ***tag density*** for all peaks within the indicated OC/CO groups at the indicated time points. Data were converted to a Z score based on the row-wise SD for each peak. ∗p \< 0.01 Mann-Whitney U test."

My question is where is this tag density information coming from.

So far for chipseq data what I have done is

* Aligned the data
* Then all required data processing post alignment
* Generated bed and bigwig files of respective samples required for comparison
* Made tag-directory using homer.
* Then peak-calling using homer,annotated peak using homer.

Now so far I am not able to figure out where is this tag-density information is .

Any help or suggestion would be really appreciated.

# Answer
> 2 votes
You can use `Homer` to retrieve tag densities, quoting from their manual:

> **Calculating ChIP-Seq Tag Densities across different experiments**
>
> annotatePeaks.pl is \[a\] useful program for cross-referencing data from multiple experiments. In order to count the number of tags from different sequencing experiments, you must first create tag directories for each of these experiments. Once created, tag counts from these directories in the vicinity of your peaks can be added by specifying "-d ...". You can specify as many tag directories as you like. Tag totals for each directory will be placed in new columns starting on column 18. For example:
>
> annotatePeaks.pl pu1peaks.txt mm8 -size 400 -d Macrophage-PU.1/ Bcell-PU.1/ \> output.txt

---
Tags: chip-seq
---


# Question
Title: Is there a way to measure cell line similarity using python?
I have a list of cell lines that should be compared against another list of cell lines as pairs. I am trying to find a way in order for me to compare two cells such as **MCF7** and **10964c** cell lines. I am using python so, can anyone provide some steps so that I can do the comparison.

I also have their ontologies in an OBO format as I need them to do the comparison. if there is a library or a certain procedure that can solve this problem will be a great help.

**Edit:**

The comparison I'm trying to do is to measure the similarity in different ways, there are many measurement methods such as the distance in the ontological hierarchy between two cells by calculating the shortest path. You can also measure the number of common ancestors. But also, there are libraries that perform complicated (and more accurate) similarity measurements like JiangConrath, Jaccard, or Resnik which is popular. The problem with theses libraries is that they are not oriented to calculate cell line ontologies and may arequire annotation corpus for some similarity measures.

# Answer
> 1 votes
At its current form, your question is not clear, for example which aspects of these cell lines are your trying to compare?

Regarding your point on the ontologies, in `Python` you can use the `pronto` library to work with ontologies in OBO format.

---
Tags: python, cell-line, ontology
---


# Question
Title: Fasterq-dump: --split-spot or -concatenate-reads?
After using files that I downloaded from the SRA with fasterq-dump, I realize I am not 100% sure that I have all the data.

I noticed in my downstream analysis that I seem to be missing the .1 and .2 numbers in the code associated with individual reads with the same spot ID. Then I noticed that my initial fastq file doesn't seem to have those either.

So I'm wondering if the -split-spot option that I used is doing what I want it do do, which is 'get all the sequence data', and how does it differ from --concatenate-reads. I am afraid I did not understand this from the program description. Please help?

# Answer
Fasterq comes from the latest version of sratools. So if you check the manual , it says the equivalence is:

> fastq-dump SRRXXXXXX --split-3 --skip-technical
>
> fasterq-dump SRRXXXXXX

In older versions of sratoolkit, if you use fastq-dump without specifying --split-3 for paired-end reads, you get the format mentioned, spotID.1 for forward, spotID.2 for reverse:

```
fastq-dump SRR6462984.sra
more SRR6462984.fastq

```

You get:

```
@SRR6462984.1 1 length=301
NATCTGCCCGTTCCACATAAACACCGTTAAACATCGATAGAGCGAAATAAAGCTGCTGAGTGATGTGATGATAAAGCTTTTCCGTTCTGACACGATCTTCCATCTCGTCTATCATTGCATCGAGAGCGTCAGAGATCGCAAGCAGTGCGGATATTGCCGTCCGTCTTGAAGCTTTGGCAGAGCCAGGCGGAATTTGTTTATCGGATAGTGTTTATGCGCAAATTTAGCATATGATTTCCGACAACTTTTGCTTCATCGGGCCGCAAAGCGTGAAAAACAGTGTCGTCACCGTTGATGTGTG
+SRR6462984.1 1 length=301
#/AAAEE/EEEA/EAAEEAEE/EE/EEA/EEEEE/AAE/EE/EAEAE//EEAEE/EEEEE<EEEEE/EE<EE/AEEEE/EAEEEEA/E<EEAEAEEEEEEAEEEE<EEAEEEEEAAAEEE/EEEEE/E//AEEEE6EA66EE<<<EAAA/AAAAAE/E6E/EA/AE/AEEEEEEAAEEEEA6EEEEEEEEE/EEE/AEEAAAEEEAE/E/EEEEEA//E//EEE/EEE/A<//AE/EE//6A///AA/E//E/E/EE//AE//AAEEE////E/A<AAE/EEAE//A//6/A</<6/</E/
@SRR6462984.2 2 length=302
NTCAGCAGGCGCCGTAACTTCAGCAAAACTGTAGCTATAACTTTTCCGAACTTTATACCTAACTTCGCTTTCCTAGCACTGTATAAAATGCCTAAAAGAAATCAGCAGGCGCCGTCACTTCAAAAAATCAGTATGGATCAATCAGACGCGGAGCTTATTAAAACAGAGCTTATTAAAAATTATGTCTTTCTTCTTCTTTCAGCAATTCGCAACGTAAAATTCAAAATTTTGGACGCCGTAATTGCGCAGGAACCAATTCAAAGAGGCGGCGGCCAATTAATTCCAACGTCTCGAAGAGTGGC
+SRR6462984.2 2 length=302

```

With the latest version, it by defaults splits the paired end sra file into forward and reverse fastq (i.e --split-3):

```
fasterq-dump SRR6462984.sra
head -2 SRR6462984_*.fastq

```

And you get

```
==> SRR6462984_1.fastq <==
@SRR6462984.1 1 length=150
NATCTGCCCGTTCCACATAAACACCGTTAAACATCGATAGAGCGAAATAAAGCTGCTGAGTGATGTGATGATAAAGCTTTTCCGTTCTGACACGATCTTCCATCTCGTCTATCATTGCATCGAGAGCGTCAGAGATCGCAAGCAGTGCGG

==> SRR6462984_2.fastq <==
@SRR6462984.1 1 length=151
ATATTGCCGTCCGTCTTGAAGCTTTGGCAGAGCCAGGCGGAATTTGTTTATCGGATAGTGTTTATGCGCAAATTTAGCATATGATTTCCGACAACTTTTGCTTCATCGGGCCGCAAAGCGTGAAAAACAGTGTCGTCACCGTTGATGTGTG

```

So if you need to get back the .1 and .2 format, use fastq-dump from the latest installation.

I would suggest sticking to the newest format, with 2 files, because fastq-dump will be deprecated in the future.
> 3 votes

---
Tags: fastq, sratoolkit
---


title: WGCNA module preservation

text:
# Question
Title: WGCNA module preservation
I am completing the online tutorials for module preservation of human brain modules in chimpanzee brains and vice versa (III). If anyone is aware of this tutorial or practiced before wondering if can help with quick question please. I completed sections 1a and 2a to upload the expression data and then proceeded to Appendix A to complete Network Analysis of expression data. Within the code; I get error: Object ConnectivityHuman or Error: Object ConnectivityChimp not found. I have searched for any Sections of tutorial missed but unable to find these? Anyone any advice on tutorial gratefully appreciated thank you

```
ConnectivityHuman = ConnectivityHuman/max(ConnectivityHuman)
    ConnectivityChimp = ConnectivityChimp/max(ConnectivityChimp)

```

# Answer
> 1 votes
Thanks for pointing this out. Looks like we forgot to include two important lines. Something along the lines of

```
ConnectivityHuman = colSums(adjHuman)-1
ConnectivityChimp = colSums(adjChimp)-1

```

This should do it; I'll update the tutorial.

---
Tags: r, wgcna
---


# Question
Title: How to calculate the number of cells in different cell clusters with Seurat v3?
I was supposed to use the code below to get the number of cells per cluster. The code works with Seurat version 2, but while using version 3 I got the error `no slot of name "ident" for this object of class "Seurat"`. Does anyone have any idea about how to deal with this?

```
# Calculate number of cells per cluster from object@ident

cell.num <- table(shox2crepos_mtmg@ident)

# Add cell number per cluster to cluster labels

ClusterLabels = paste("Cluster", names(cell.num), paste0("(n = ", cell.num, ")"))

# Order legend labels in plot in the same order as 'ClusterLabels'

ClusterBreaks = names(cell.num)

# Plot tSNE with new legend labels for clusters

TSNEPlot(object = shox2crepos_mtmg, do.return = T) +
  scale_colour_discrete(breaks = ClusterBreaks,
                        labels = ClusterLabels) +
  labs(x = "t-SNE 1",
       y = "t-SNE 2")

```

# Answer
The "identities" can be accessed with the `Idents()` function in Seurat v3.

You should change `cell.num <- table(shox2crepos_mtmg@ident)` to

`cell.num <- table(Idents(shox2crepos_mtmg))`
> 2 votes

---
Tags: scrnaseq
---


# Question
Title: Summing two columns where they share interval overlap
I have two bed files where each column has a numeric score. I want to sum these numbers where there is overlap between the two bed files. How do I do that? Also, I am only interested in the score in specific genes.

Here is my setup:

```
import numpy as np
np.random.seed(0)
import pyranges as pr
import pandas as pd

entries = int(1e6)
cpg1 = pr.random(entries)
cpg2 = pr.random(entries)

cpg1.CPG = np.random.randint(100, size=entries)
cpg2.CPG = np.random.randint(100, size=entries)

regulatory_elements = pr.random(int(1e5), length=int(1e4))
regulatory_elements.Gene = np.arange(len(regulatory_elements), dtype=int)

```

So my data looks like:

```
cpg1
# +--------------+-----------+-----------+--------------+-----------+
# | Chromosome   | Start     | End       | Strand       | CPG       |
# | (category)   | (int32)   | (int32)   | (category)   | (int64)   |
# |--------------+-----------+-----------+--------------+-----------|
# | chr1         | 8830650   | 8830750   | +            | 72        |
# | chr1         | 9564361   | 9564461   | +            | 57        |
# | chr1         | 44977425  | 44977525  | +            | 99        |
# | chr1         | 239741543 | 239741643 | +            | 54        |
# | ...          | ...       | ...       | ...          | ...       |
# | chrY         | 29437476  | 29437576  | -            | 44        |
# | chrY         | 49995298  | 49995398  | -            | 43        |
# | chrY         | 50840129  | 50840229  | -            | 10        |
# | chrY         | 38069647  | 38069747  | -            | 89        |
# +--------------+-----------+-----------+--------------+-----------+
# Stranded PyRanges object has 1,000,000 rows and 5 columns from 25 chromosomes.
# For printing, the PyRanges was sorted on Chromosome and Strand.

regulatory_elements
# +--------------+-----------+-----------+--------------+-----------+
# | Chromosome   | Start     | End       | Strand       | Gene      |
# | (category)   | (int32)   | (int32)   | (category)   | (int64)   |
# |--------------+-----------+-----------+--------------+-----------|
# | chr1         | 217921634 | 217931634 | +            | 0         |
# | chr1         | 166226804 | 166236804 | +            | 1         |
# | chr1         | 170688210 | 170698210 | +            | 2         |
# | chr1         | 57958563  | 57968563  | +            | 3         |
# | ...          | ...       | ...       | ...          | ...       |
# | chrY         | 35870043  | 35880043  | -            | 99996     |
# | chrY         | 35057634  | 35067634  | -            | 99997     |
# | chrY         | 52718975  | 52728975  | -            | 99998     |
# | chrY         | 19196507  | 19206507  | -            | 99999     |
# +--------------+-----------+-----------+--------------+-----------+
# Stranded PyRanges object has 100,000 rows and 5 columns from 25 chromosomes.
# For printing, the PyRanges was sorted on Chromosome and Strand.

```

(This is a q from privateish correspondence).

# Answer
If you want to sum two columns based on overlaps, converting them to RLEs is a good idea.

```
# get the cpg scores only where there is overlap with the genes you are interested in
cr1 = cpg1.join(regulatory_elements).new_position("intersection").drop(like="_(a|b)")
cr2 = cpg2.join(regulatory_elements).new_position("intersection").drop(like="_(a|b)")

# convert these to run length encodings, using the score from the CPG column
rle1 = cr1.to_rle("CPG", rpm=False)
rle2 = cr2.to_rle("CPG", rpm=False)

rle_sum = rle1 + rle2

# convert the rles back to ranges
# (I will give the to_ranges function an option dtype eventually)
range_sum = rle_sum.to_ranges().apply(lambda df: df.astype({"Start": np.int32, "End": np.int32}))

# add the gene names to the scores
gr = range_sum.join(regulatory_elements) # .drop(like="_b")

```

Now the result is:

```
print(gr)
# +--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------+
# | Chromosome   | Start     | End       | Score       | Strand     | Start_b   | End_b     | Strand_b     | Gene      |
# | (object)     | (int32)   | (int32)   | (float64)   | (object)   | (int32)   | (int32)   | (category)   | (int64)   |
# |--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------|
# | chr1         | 15623     | 15723     | 77.0        | +          | 9886      | 19886     | +            | 3714      |
# | chr1         | 16034     | 16134     | 84.0        | +          | 9886      | 19886     | +            | 3714      |
# | chr1         | 36588     | 36688     | 2.0         | +          | 36479     | 46479     | +            | 2871      |
# | chr1         | 37271     | 37371     | 30.0        | +          | 36479     | 46479     | +            | 2871      |
# | ...          | ...       | ...       | ...         | ...        | ...       | ...       | ...          | ...       |
# | chrY         | 59312851  | 59312951  | 28.0        | -          | 59308107  | 59318107  | +            | 98599     |
# | chrY         | 59332843  | 59332943  | 82.0        | -          | 59332664  | 59342664  | -            | 99215     |
# | chrY         | 59337592  | 59337692  | 41.0        | -          | 59332664  | 59342664  | -            | 99215     |
# | chrY         | 59341719  | 59341819  | 15.0        | -          | 59332664  | 59342664  | -            | 99215     |
# +--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------+
# Stranded PyRanges object has 668,317 rows and 9 columns from 25 chromosomes.
# For printing, the PyRanges was sorted on Chromosome and Strand.

print(gr[gr.Score > 100])
# +--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------+
# | Chromosome   | Start     | End       | Score       | Strand     | Start_b   | End_b     | Strand_b     | Gene      |
# | (object)     | (int32)   | (int32)   | (float64)   | (object)   | (int32)   | (int32)   | (category)   | (int64)   |
# |--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------|
# | chr1         | 96874     | 96974     | 116.0       | +          | 87668     | 97668     | -            | 7391      |
# | chr1         | 96874     | 96974     | 116.0       | +          | 89794     | 99794     | -            | 6357      |
# | chr1         | 749686    | 749697    | 180.0       | +          | 744041    | 754041    | +            | 1150      |
# | chr1         | 749686    | 749697    | 180.0       | +          | 749277    | 759277    | -            | 6395      |
# | ...          | ...       | ...       | ...         | ...        | ...       | ...       | ...          | ...       |
# | chrY         | 59060667  | 59060767  | 196.0       | -          | 59052942  | 59062942  | -            | 99957     |
# | chrY         | 59060667  | 59060767  | 196.0       | -          | 59059904  | 59069904  | +            | 98401     |
# | chrY         | 59063570  | 59063576  | 105.0       | -          | 59059904  | 59069904  | +            | 98401     |
# | chrY         | 59269010  | 59269085  | 152.0       | -          | 59264704  | 59274704  | -            | 99499     |
# +--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------+
# Stranded PyRanges object has 106,643 rows and 9 columns from 25 chromosomes.
# For printing, the PyRanges was sorted on Chromosome and Strand.

```
> 0 votes

---
Tags: python, pybedtools, pyranges
---


# Question
Title: Range overlap python error with genomic regions
I have two files

```
s3.txt :
1   10  20
1   5   20
2   20  30
2   25  30
1   10  50
2   20  60
1   14  17

s4.txt:
1   10  20
2   20  30

```

I am trying to match col0 of both the files and get rows that fall between range(inclusive of themselves) 10-20 and 20-30 as seen in s4 file. file s4 has co ordinates which can be used as reference range (chrom start and end) and s3 has list of co ordinates from an experimental condition what I am trying to achieve is to which co-ordinates from my file s3 fall on or between my reference co ordinates in s4.

code so far:

```
containing_ranges = []
with open('s4.txt', 'r') as f:
    for line in f:
        fields = line.strip().split('\t')
        containing_ranges.append(fields)

tested_ranges = []
with open('s3.txt', 'r') as f:
    for line in f:
        fields = line.strip().split('\t')
        tested_ranges.append(fields)

for c_range in containing_ranges:
    for t_range in tested_ranges:
        tst = int(t_range[1])
        ten = int(t_range[2])
        cst = int(c_range[1])
        cen = int(c_range[2])
        if  c_range[0] == t_range[0]:
            included = cst >= tst and cen <= ten
            if included == True:
               print t_range

```

Output with missing row(1 14 17) :

```
['1', '10', '20']
['1', '5', '20']
['1', '10', '50']
['2', '20', '30']
['2', '20', '60']

```

Desired output:

```
1   10  20
2   20  30
2   25  30
1   14  17

```

Not sure if my logic is wrong and why does it miss 14-17 as it falls between 10-20

```
[EDIT] using pybedtools

>>> print(s4.intersect(s3, wb=True))
1   10  20  1   10  20
1   10  20  1   5   20
1   10  20  1   10  50
1   14  17  1   14  17
2   20  30  2   20  30
2   25  30  2   25  30
2   20  30  2   20  60

>>> print(s4.intersect(s3, wa=True, wb=True, F=1))
1   10  20  1   10  20
1   10  20  1   14  17
2   20  30  2   20  30
2   20  30  2   25  30

using bedops
bin$ less answer.bed
1       5       20
1       10      20
1       10      50
1       14      17
2       20      30
2       20      60
2       25      30

using @bli code(on python2.7)
('1', 10, 20)
('1', 14, 17)
('2', 20, 30)
('2', 25, 30)
 why can I not see the interval 1 5 20

```

# Answer
> 1 votes
pyranges answer:

```
# pip install pyranges
# or
# conda install -c bioconda pyranges

import pyranges as pr
s3 = pr.read_bed("s3.bed")
s4 = pr.read_bed("s4.bed")
s3.intersect(s4, how="containment")

```

---

Answer with setup (Edit):

```
import pandas as pd
from io import StringIO
import pyranges as pr

c1 = """1   10  20
1   5   20
2   20  30
2   25  30
1   10  50
2   20  60
1   14  17"""

c2 = """1   10  20
2   20  30"""

columns = "Chromosome Start End".split()

df1 = pd.read_table(StringIO(c1), sep="\s+", header=None, names=columns)
df2 = pd.read_table(StringIO(c2), sep="\s+", header=None, names=columns)

gr1 = pr.PyRanges(df1)
gr2 = pr.PyRanges(df2)

print(gr1.intersect(gr2, how="containment"))

```

Result:

```
+--------------+-----------+-----------+
|   Chromosome |     Start |       End |
|   (category) |   (int32) |   (int32) |
|--------------+-----------+-----------|
|            1 |        10 |        20 |
|            1 |        14 |        17 |
|            2 |        20 |        30 |
|            2 |        25 |        30 |
+--------------+-----------+-----------+
Unstranded PyRanges object has 4 rows and 3 columns from 2 chromosomes.
For printing, the PyRanges was sorted on Chromosome.

```

# Answer
> 5 votes
You're reinventing `bedtools intersect` (or bedops), for which there's already a convenient python module:

```
from pybedtools import BedTool

s3 = BedTool('s3.bed')
s4 = BedTool('s4.bed')

print(s4.intersect(s3, wa=True, wb=True, F=1))

```

The `wb=True` is equivalent to `-wb` with `bedtools intersect` on the command line. Similarly, `F=1` is the same as `-F 1`.

# Answer
> 3 votes
### Regarding your python code

If you want the **experimental ranges that are entirely contained in one of the reference ranges**, you need to have the coordinates in the following order:

```
cst <= tst < ten <= cen

```

If what you want are the **experimental ranges that overlap one of the reference ranges**, you need to have either the start or the end of the experimental range fall within the reference range:

```
(cst <= tst < cen) or (cst < ten <= cen)

```

Your code is equivalent to neither possibilities:

```
cst >= tst and cen <= ten

```

This is equivalent to:

```
tst <= cst and cen <= ten

```

(Or `tst <= cst < cen <= ten`, since we know that `cst < cen`, by definition of a bed interval).

With this rewriting, we can more easily see that you are actually selecting the **experimental ranges that contain a reference range**.

Here is some (python3) code that gives you the results in the other two situations:

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

ref_intervals = []
with open("s4.txt", "r") as f:
    for line in f:
        (chr, start, end) = line.strip().split("\t")
        ref_intervals.append((chr, int(start), int(end)))

exp_intervals = []
with open("s3.txt", "r") as f:
    for line in f:
        (chr, start, end) = line.strip().split("\t")
        exp_intervals.append((chr, int(start), int(end)))

contained = []
overlapping = []
for (r_chr, r_start, r_end) in ref_intervals:
    for (e_chr, e_start, e_end) in exp_intervals:
        if e_chr == r_chr:
            if r_start <= e_start < r_end or r_start < e_end <= r_end:
                overlapping.append((e_chr, e_start, e_end))
            if r_start <= e_start < e_end <= r_end:
                contained.append((e_chr, e_start, e_end))

print("overlapping")
for (chr, start, end) in overlapping:
    print(chr, start, end, sep="\t")

print("contained")
for (chr, start, end) in contained:
    print(chr, start, end, sep="\t")

```

If I run it, I obtain the following results:

```
$ ./overlap.py
overlapping
1   10  20
1   5   20
1   10  50
1   14  17
2   20  30
2   25  30
2   20  60
contained
1   10  20
1   14  17
2   20  30
2   25  30

```

It is probably an good exercise to program this, but as the other answers point out, **there are efficient tools that would be a better solution in a professional setting**.

# Answer
> 1 votes
You could use BEDOPS, instead:

```
$ sort-bed s3.txt > s3.bed
$ sort-bed s4.txt > s4.bed
$ bedops --element-of 1 s3.bed s4.bed > answer.bed

```

If you need to run it from within Python, you could use `subprocess.check_output()`:

```
import subprocess
...
try:
    result = subprocess.check_output("bedops --element-of 1 %s %s > %s" % (set_a_fn, set_b_fn, answer_fn), shell=True)
except subprocess.CalledProcessError as err:
    raise SystemExit("Could not run bedops\n")
# do stuff with 'result'

```

---
Tags: bed, genomics, python, homework
---


# Question
Title: ATAC-seq macs2 peak splitting in sliding windows
*This question has also been asked on Biostars*

I used macs2 to call peaks for atac-seq data. now my goal is to split the peaks into 50 bp windows with 25 bp steps and then calculate the Tn5 integration frequency in each window.

How should I proceed with that?

# Answer
```
import numpy as np
np.random.seed(0)
import pyranges as pr

gr = pr.random()
gr.Score = np.random.randint(100, size=len(gr))

gr = gr.slack(25) # make data wider for this example

print(gr)

t1 = gr.tile(50)

def increase_by_25(df):
    df = df.copy()
    df.Start += 25
    df.End += 25
    return df

t2 = t1.apply(increase_by_25)

tiled = pr.concat([t1, t2]).sort()

print(tiled)
# +--------------+-----------+-----------+--------------+-----------+
# | Chromosome   | Start     | End       | Strand       | Score     |
# | (category)   | (int32)   | (int32)   | (category)   | (int64)   |
# |--------------+-----------+-----------+--------------+-----------|
# | chr1         | 5205300   | 5205350   | +            | 33        |
# | chr1         | 5205325   | 5205375   | +            | 33        |
# | chr1         | 5205350   | 5205400   | +            | 33        |
# | chr1         | 5205375   | 5205425   | +            | 33        |
# | ...          | ...       | ...       | ...          | ...       |
# | chrY         | 41326450  | 41326500  | -            | 3         |
# | chrY         | 41326475  | 41326525  | -            | 3         |
# | chrY         | 41326500  | 41326550  | -            | 3         |
# | chrY         | 41326525  | 41326575  | -            | 3         |
# +--------------+-----------+-----------+--------------+-----------+
# Stranded PyRanges object has 7,964 rows and 5 columns from 24 chromosomes.
# For printing, the PyRanges was sorted on Chromosome and Strand.

```

This is as far as I can get without example data and a clearer explanation of what you need.
> 1 votes

---
Tags: atac-seq, macs2
---


# Question
Title: What type of information does an .sf file contain?
What is an .sf file. What information does it contain?

# Answer
You should provide more context, otherwise one can only guess what is it you are looking for! In case you are asking about `Salmon`'s output, below is from their documentation:

> Salmon’s main output is its quantification file. This file is a plain-text, tab-separated file with a single header line (which names all of the columns). This file is named quant.sf and appears at the top-level of Salmon’s output directory. The columns appear in the following order:

```
Name   Length   EffectiveLength   TPM   NumReads

```
> 0 votes

---
Tags: file-formats
---


# Question
Title: Are there genomics databases for neurological and psychiatric disorders?
Are there databases that relate genetic disorders with neurological disorders?

I am thinking of something like COSMICS for cancer? AFAIK cancer usually develops due to a combination of 5-10 mutations.

Also if you can comment how genomics/bioinformatics in neurological contexts is compared to the how bioinformaticians study cancer?

Related: https://www.nimhgenetics.org/resources/psychencode

# Answer
> 3 votes
Try with PsyGeNet:

> PsyGeNET (Psychiatric disorders Gene association NETwork) is a resource for the exploratory analysis of psychiatric diseases and their associated genes

It is based on literature mining and curated by experts.

Here is the publication if you want more details:

PsyGeNET: a knowledge platform on psychiatric disorders and their genes

and they even have a Bioconductor package for programatic queries:

psygenet2r: a R/Bioconductor package for the analysis of psychiatric disease genes

---
Tags: database, genomics
---


title: Loading local FASTA file in igv.js

text:
# Question
Title: Loading local FASTA file in igv.js
I'm trying to use igv.js for a very simple visualization of some genomic alignments, but I can't get past the first step. I have a FASTA file that I want to use as the reference, but igv.js consistently gives an error when I try to load that file. The file is a valid FASTA (have tested others from other sources to be sure - same error), and I can specify a remote https:// file, but not a local file.

Here's the HTML page I'm working on:

```
<!DOCTYPE html>
<html lang="en">

<head>
    <meta charset="UTF-8">
    <meta name="viewport" content="width=device-width, initial-scale=1, maximum-scale=1, user-scalable=no">

    <title>Local FASTA Test</title>

    <!-- IGV JS-->
    <script src="https://igv.org/web/release/2.1.0/dist/igv.min.js"></script>

</head>

<body>

    <h1>Local File Test</h1>

    <div id="igv-div" style="padding-top: 10px;padding-bottom: 10px; border:1px solid lightgray"></div>

    <script type="text/javascript">
        var igvDiv = document.getElementById("igv-div");
        var options = {
            reference: {
                fastaURL: "test.fasta"
            }
        };

        igv.createBrowser(igvDiv, options)
            .then(function(browser) {
                console.log("Created IGV browser");
            })
    </script>

</body>

</html>

```

and here's the FASTA file

```
>HSBGPG Human gene for bone gla protein (BGP)
GGCAGATTCCCCCTAGACCCGCCCGCACCATGGTCAGGCATGCCCCTCCTCATCGCTGGGCACAGCCCAGAGGGT
ATAAACAGTGCTGGAGGCTGGCGGGGCAGGCCAGCTGAGTCCTGAGCAGCAGCCCAGCGCAGCCACCGAGACACC
ATGAGAGCCCTCACACTCCTCGCCCTATTGGCCCTGGCCGCACTTTGCATCGCTGGCCAGGCAGGTGAGTGCCCC
CACCTCCCCTCAGGCCGCATTGCAGTGGGGGCTGAGAGGAGGAAGCACCATGGCCCACCTCTTCTCACCCCTTTG
GCTGGCAGTCCCTTTGCAGTCTAACCACCTTGTTGCAGGCTCAATCCATTTGCCCCAGCTCTGCCCTTGCAGAGG
GAGAGGAGGGAAGAGCAAGCTGCCCGAGACGCAGGGGAAGGAGGATGAGGGCCCTGGGGATGAGCTGGGGTGAAC
CAGGCTCCCTTTCCTTTGCAGGTGCGAAGCCCAGCGGTGCAGAGTCCAGCAAAGGTGCAGGTATGAGGATGGACC
TGATGGGTTCCTGGACCCTCCCCTCTCACCCTGGTCCCTCAGTCTCATTCCCCCACTCCTGCCACCTCCTGTCTG
GCCATCAGGAAGGCCAGCCTGCTCCCCACCTGATCCTCCCAAACCCAGAGCCACCTGATGCCTGCCCCTCTGCTC
CACAGCCTTTGTGTCCAAGCAGGAGGGCAGCGAGGTAGTGAAGAGACCCAGGCGCTACCTGTATCAATGGCTGGG
GTGAGAGAAAAGGCAGAGCTGGGCCAAGGCCCTGCCTCTCCGGGATGGTCTGTGGGGGAGCTGCAGCAGGGAGTG
GCCTCTCTGGGTTGTGGTGGGGGTACAGGCAGCCTGCCCTGGTGGGCACCCTGGAGCCCCATGTGTAGGGAGAGG
AGGGATGGGCATTTTGCACGGGGGCTGATGCCACCACGTCGGGTGTCTCAGAGCCCCAGTCCCCTACCCGGATCC
CCTGGAGCCCAGGAGGGAGGTGTGTGAGCTCAATCCGGACTGTGACGAGTTGGCTGACCACATCGGCTTTCAGGA
GGCCTATCGGCGCTTCTACGGCCCGGTCTAGGGTGTCGCTCTGCTGGCCTGGCCGGCAACCCCAGTTCTGCTCCT
CTCCAGGCACCCTTCTTTCCTCTTCCCCTTGCCCTTGCCCTGACCTCCCAGCCCTATGGATGTGGGGTCCCCATC
ATCCCAGCTGCTCCCAAATAAACTCCAGAAG

```

# Answer
> 2 votes
It looks as though the script is looking for a URL (fastaURL) so you could try specifying the file as file://directory/test.fasta

By default though, javascript cannot read local files from disk for security reasons. Are you sure it is even possible to specify a local file in your script in this way? If yes, then the problem is likely that you need to specify the path correctly.

# Answer
> 1 votes
just stumbled on this. I'm the developer of igv.js. Please feel free to open issues and questions on our git forum, I don't often have time to browse other forums.

Browsers will not allow loading of local files by file path in javascript, if they did it would be a gigantic security hole. This restriction includes "file://" urls. Local files can be loaded by user action (clicking a button), that is supported in igv.js but there is no button for it as igv.js is designed as a component to embed in pages, not an application by itself. For an example see igv.org/app (https://github.com/igvteam/igv-webapp), a minimal application wrapping igv.js. You'll find "local file" options for genomes and tracks there.

---
Tags: visualization, igv
---


# Question
Title: Is re-normalization of RNAseq data recommended for analysis of gene subsets?
I downloaded an RNAseq dataset from TCGA database in 3 formats: 1) HTSeq counts; 2) FPKM; 3) FPQM-upper quartile normalized.

The complete dataset contains ~60,000 genes. All of my analysis will focus on a subset of ~2500 genes.

What confused me is the following statement on the TCGA documentation page: "RNA-Seq expression level read counts produced by HT-Seq are normalized using two similar methods: FPKM and FPKM-UQ. Normalized values should be used only within the context of the entire gene set. Users are encouraged to normalize raw read count values if a subset of genes is investigated."

If I normalize the 2500 gene counts in isolation of the rest, I am guessing the normalization factor will be more unstable between samples than using the complete gene set. It means that if one patient has higher expression of the subset than another, normalization within subset will obscure that difference.

However, the normalization factor (either total mapped reads or upper quartile) becomes more stable the larger the number of genes considered. Less variation will be seen on a global scale because of biological differences, whilst technical differences will be normalized.

In short, do you normalize the gene counts de novo if you analyse only a subset of genes? Why?

# Answer
FPKM, RPKM are ways of normalising for both sequencing depth and gene length. You first calculate a scaling factor (`SF = #read/1Mil`),you then divide the genes counts for the `SF`(`RPM = # counts/SF`), and finally you divide by the length of the gene (`RPM/length(gene)`). Each sample has a different `SF` meaning that you should not compare values from different samples.

To fix this, you have to normalise across samples, by `upper quartile normalise` using a subset of genes. In this case, the `SF` is calculated based on the upper quartile (70%) of this shared and immutable subset, and it is therefore the same for all the samples.

If you subset your genes, and you are using the normalised values for downstream analyses, then you need to re-scale your values.

I personally try to avoid `FPKM`,`RPKM`,`TPM` (if I can) and always start from raw counts.

I would suggest to have a look at the documentation of packages like `DESeq2` and `edgeR` for detailed and better explanations.
> 2 votes

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


title: CIBERSORT runtime error

text:
# Question
Title: CIBERSORT runtime error
I was running CIBERSORT web and saw this error:

```
Exception in thread "Thread-0" java.lang.ArrayIndexOutOfBoundsException: 0
at launch.QN.doQN(QN.java:56)
at launch.Cibersort.Run_CIBERSORT(Cibersort.java:611)
at launch.Cibersort.run_all(Cibersort.java:179)
at launch.Cibersort$runCIBERSORT.run(Cibersort.java:1585)
at java.lang.Thread.run(Thread.java:745)
311

```

Does anyone know why?

# Answer
> 0 votes
OK I found the bug. CIBERSORT seems to require tsv files instead of csv files. Plus it doesn't seem to like empty columns

---
Tags: webservice, cibersort
---


title: What do these files / annotations mean?

text:
# Question
Title: What do these files / annotations mean?
I have no experience in Bioinformatics and I need to understand what the annotations given here mean (I am including the first few lines, please see the link for more):

```
contig_id                         feature_id              type       location                                 start   stop    strand
gi|1045318032|gb|MAZB01000333.1|  fig|252393.25.peg.1     peg        gi|1045318032|gb|MAZB01000333.1|_511_2   511     2       -
gi|1045318033|gb|MAZB01000332.1|  fig|252393.25.repeat.1  repeat     gi|1045318033|gb|MAZB01000332.1|_1_128   1       128     +
gi|1045318033|gb|MAZB01000332.1|  fig|252393.25.repeat.2  repeat     gi|1045318033|gb|MAZB01000332.1|_318_532 318     532     +
gi|1045318033|gb|MAZB01000332.1|  fig|252393.25.peg.2     peg        gi|1045318033|gb|MAZB01000332.1|_530_321 530     321     -
gi|1045318034|gb|MAZB01000331.1|  fig|252393.25.repeat.3  repeat     gi|1045318034|gb|MAZB01000331.1|_1_534   1       534     +
gi|1045318034|gb|MAZB01000331.1|  fig|252393.25.peg.3     peg        gi|1045318034|gb|MAZB01000331.1|_3_410   3       410     +
gi|1045318035|gb|MAZB01000330.1|  fig|252393.25.repeat.4  repeat     gi|1045318035|gb|MAZB01000330.1|_1_128   1       128     +
gi|1045318035|gb|MAZB01000330.1|  fig|252393.25.peg.4     peg        gi|1045318035|gb|MAZB01000330.1|_3_539   3       539     +
gi|1045318035|gb|MAZB01000330.1|  fig|252393.25.repeat.5  repeat     gi|1045318035|gb|MAZB01000330.1|_413_539 413     539     +

```

Could someone help me to understand what these files and all annotations mean?

# Answer
You are looking at the supplementary data of a paper. That seems to have given you a list of features, and some information about those features. Specifically, you seem to have a list of two types of element:

1. "peg". Based on the information there, I assume "peg" stands for "protein encoding gene". Note that all peg lines have a peptide sequence (not seen in your question but present in the original link) and their "function" column is always some sort of protein function (e.g. "hypothetical protein" or "Transposase").
2. "repeat". These are presumably repetitive elements.

For each line, you are given an identifier for the sequence in which this element can be found. For example, your first line is on `gi|1045318032|gb|MAZB01000333.1|`. That is a GI (`gi|1045318032`) and a GenBank identifier (`gb|MAZB01000333.1`). If you go to the NCBI website, and search for `1045318032`, you will find that it is the Candidatus Erwinia dacicola isolate Erw\_SC contig\_333, whole genome shotgun sequence.

So, your first line is telling you that the authors have identified a protein encoding gene on sequence `gi|1045318032` and the gene starts at position 511 of that sequence and ends at position 2. This might sound odd, but that's because the gene is on the negative (`-`) strand and positions are calculated with respect to the positive strand. So, if you were to count on the negative, the gene would start at 2 and end at 511.

The same goes for all the other lines you have there. Just look up the GI number and you can know what sequences they are referring to.

Finally, the second column, the "feature\_id" looks like something specific to this figshare service that is hosting the data. I dont' know what it is, but I would guess it's an identifier for this particular feature. I don't think it has any relevance outside http://figshare.com.
> 5 votes

# Answer
gi looks like a Genbank Identifier. You can search for the 10 digit number following gi| at the Genbank website.
> 1 votes

---
Tags: file-formats, edger, identifiers
---


title: CIBERSORT runtime error eval failed

text:
# Question
Title: CIBERSORT runtime error eval failed
I was running CIBERSORT but caught this error right after it started permutation:

```
org.rosuda.REngine.Rserve.RserveException: eval failed, request status: error code: 127
88

```

# Answer
Ironically, I found that this problem was caused by the fact that the tables are too large! They really should have designed the correct error messages!
> 1 votes

---
Tags: rna-seq, cibersort
---


# Question
Title: Parsing MSigDB Supplementary Collections
MSigDB recently added Supplementary Collections. The formatting is similar to the original MSigDB gene sets, but there is an additional `metadata.txt` file (link).

It looks like (condensed for easier reading):

```
STANDARD_NAME   Fan_Embryonic_CTX_Big_Groups_Cajal_Retzius
ORGANISM    Homo sapiens
ORGAN_SYSTEM    Central Nervous System
PMID    29867213
PUBLICATION_TITLE   Spatial transcriptomic survey of human embryonic cerebral cortex by single-cell RNA-seq analysis
AUTHORS Fan X,Dong J,Zhong S,Wei Y,Wu Q,Yan L,Yong J,Sun L,Wang X,Zhao Y,Wang W,Yan J,Wang X,Qiao J,Tang F,
GEOID   GSE103723
...
STANDARD_NAME   Fan_Embryonic_CTX_Big_Groups_Brain_Endothelial
ORGANISM    Homo sapiens
ORGAN_SYSTEM    Central Nervous System
PMID    29867213
PUBLICATION_TITLE   Spatial transcriptomic survey of human embryonic cerebral cortex by single-cell RNA-seq analysis
AUTHORS Fan X,Dong J,Zhong S,Wei Y,Wu Q,Yan L,Yong J,Sun L,Wang X,Zhao Y,Wang W,Yan J,Wang X,Qiao J,Tang F,
GEOID   GSE103723
...
STANDARD_NAME   Fan_Embryonic_CTX_Big_Groups_Excitatory_Neuron
ORGANISM    Homo sapiens
ORGAN_SYSTEM    Central Nervous System
PMID    29867213
PUBLICATION_TITLE   Spatial transcriptomic survey of human embryonic cerebral cortex by single-cell RNA-seq analysis
AUTHORS Fan X,Dong J,Zhong S,Wei Y,Wu Q,Yan L,Yong J,Sun L,Wang X,Zhao Y,Wang W,Yan J,Wang X,Qiao J,Tang F,
GEOID   GSE103723
...

```

I could not find any description of the formatting. Is this a standard format? Is there a parser for it (preferably in R)? It looks like I can just transpose each record to get it in a table format, but I don't know if I can trust all records to be consistent.

# Answer
> 2 votes
If you want to separate each set, you can do something like this

```
dat <- readLines("~/Desktop/scsig.v1.0.metadata.txt")
dat <- gsub("\t"," ",dat)

# get the indeces of each type
standardName <- grep("STANDARD_NAME",dat)
organism <- grep("ORGANISM",dat)
organSystem <- grep("ORGAN_SYSTEM",dat)
pmid <- grep("PMID",dat)

# return the values
standardName <- gsub(pattern = "STANDARD_NAME",replacement = "",dat[standardName])
organism <- gsub(pattern = "ORGANISM",replacement = "",dat[organism])
organSystem <- gsub(pattern = "ORGAN_SYSTEM",replacement = "",dat[organSystem])
pmid <- gsub(pattern = "PMID",replacement = "",dat[pmid])

df <- data.frame(standardName,organism,organSystem,pmid,stringsAsFactors = FALSE)

 head(df)
                                     standardName      organism             organSystem      pmid
1      Fan_Embryonic_CTX_Big_Groups_Cajal_Retzius  Homo sapiens  Central Nervous System  29867213
2  Fan_Embryonic_CTX_Big_Groups_Brain_Endothelial  Homo sapiens  Central Nervous System  29867213
3  Fan_Embryonic_CTX_Big_Groups_Excitatory_Neuron  Homo sapiens  Central Nervous System  29867213
4              Fan_Embryonic_CTX_Big_Groups_Glial  Homo sapiens  Central Nervous System  29867213
5       Fan_Embryonic_CTX_Big_Groups_Brain_Immune  Homo sapiens  Central Nervous System  29867213
6          Fan_Embryonic_CTX_Big_Groups_Microglia  Homo sapiens  Central Nervous System  29867213

```

This method is very verbose, but it's a start and it can be easily improved with a loop.

```
 #---- EDIT

```

Not knowing how the format works, there may be two issues:

1) gene sets with missing values only

```
STANDARD_NAME Aizarani_Liver_C33_Stellate_cells_2
ORGANISM
ORGAN_SYSTEM
PMID 31292543

```

2) gene sets with missing value and type

```
STANDARD_NAME Aizarani_Liver_C33_Stellate_cells_2
PMID 31292543

```

`Case 1` is not an issue, as the code above will simply add an empty cell and go to the next line. For `Case 2` we can use `unique` to get a list of all define types, create an empty `data.frame` with those `types` and have a loop (or `sapply`) to add the values in the right columns.

We can check the nature of `metadata.txt`:

```
 types <- gsub("^(.*)\t.*$","\\1",dat)
 sapply(unique(types), function(t) length(grep(t,types)))

     STANDARD_NAME                ORGANISM            ORGAN_SYSTEM                    PMID       PUBLICATION_TITLE
                257                     257                     257                     257                     257
            AUTHORS                   GEOID            EXACT_SOURCE    EXTERNAL_DETAILS_URL                    CHIP
                257                     257                     257                     257                     257
      CATEGORY_CODE             CONTRIBUTOR         CONTRIBUTOR_ORG       DESCRIPTION_BRIEF RAW_PUBLICATION_MEMBERS
                257                     514                     257                     257                     257

```

Your are lucky and you don't have to worry about shifting because of missing values. There is a fixable issue with `CONTRIBUTORS`, but it is related to the `regex`, not the file format.

---
Tags: r, file-formats, format-conversion, pathway
---


# Question
Title: Customizing colors in geom_arc_bar (ggforce)
I am trying to plot scatter plot of pie charts using ggplot2 and ggforce. Following is my code

```
library(data.table)
library(ggplot2)
library(ggforce)
library(dplyr)

df <- fread("data_frame_pie_plot.csv", header = TRUE,
        data.table = FALSE, stringsAsFactors = FALSE)

req_data <- apply(df, 1, function(x){
  geom_arc_bar(aes(x0 = x0, y0 = y0, r0 = 0, r = 0.45,
               amount = amount, fill = discrete_fill), stat = "pie")})

ggplot(df) + req_data + theme_no_axes()

```

plot\_obtained Link\_to\_data Data\_frame.
Can I fill the colours manually (discrete(based on discrete\_fill column) and continous (based on conti\_fill column))? I have tried the following codes.

```
##================================================================================##
# Code 1
my_colors <- c("red","white","blue")
ggplot(df, fill = sample_mut_colour) + req_data +
  theme_no_axes() + scale_color_manual(values=my_colors)

##================================================================================##
# Code 2
my_colors_1 <- c("red", "white")
my_colors_2 <- c("blue", "white")

req_data <- apply(df, 1, function(x){
  geom_arc_bar(aes(x0 = x0, y0 = y0, r0 = 0, r = 0.45,
               amount = amount,
               fill = discrete_fill,
               color = case_when(discrete_fill == 1 ~ my_colors_1,
                                 discrete_fill == 2 ~ my_colors_2)), stat = "pie")})

ggplot(df) + req_data + theme_no_axes()

##================================================================================##
# Code 3
req_data <- apply(df, 1, function(x){
  geom_arc_bar(aes(x0 = x0, y0 = y0, r0 = 0, r = 0.45,
               amount = amount,
               fill = case_when(discrete_fill == 1 ~ my_colors_1,
                                 discrete_fill == 2 ~ my_colors_2)), stat = "pie")})

ggplot(df) + req_data + theme_no_axes()

##================================================================================##
# Code 4
req_data <- apply(df, 1, function(x){
  geom_arc_bar(aes(x0 = x0, y0 = y0, r0 = 0, r = 0.45,
               amount = amount,
               fill = discrete_fill),
           palette = case_when(discrete_fill == 1 ~ "red",
                               discrete_fill == 2 ~ "black",
                               TRUE ~ "white"), stat = "pie")})

ggplot(df) + req_data + theme_no_axes()

##=================================================================================##

```

But they do not seem to work. Any idea what is wrong in the code?

# Answer
> 2 votes
For a discrete scale your Code 1 example was nearly there. Just use `scale_fill_manual` instead of `scale_color_manual`. For more control over how the discrete values map to color values you can name the vector, too.

```
my_colors <- c("red","white","blue")
ggplot(df) + req_data +
    theme_no_axes() + scale_fill_manual(values=my_colors)
# or, if you want consistent output independent of the order of
# df$discrete_fill:
my_colors <- c("1" = "red", "2" = "white", "dummy" = "blue")
ggplot(df) + req_data +
    theme_no_axes() + scale_fill_manual(values=my_colors)

```

You can also just use a single call to geom\_arc\_bar for the whole data frame (rather than a list of call results from each row) and you'll get the same effect.

```
ggplot(df) +
  geom_arc_bar(aes(x0 = x0, y0 = y0, r0 = 0, r = 0.45,
               amount = amount, fill = discrete_fill), stat = "pie") +
    theme_no_axes() + scale_fill_manual(values=my_colors)

```

---
Tags: r, ggplot2
---


# Question
Title: bed file with N regions of GRCh38 reference?
How do I obtain a bed file with the list of N-regions of GRCh38 reference? This is, the regions where the sequence is a stretch of Ns.

# Answer
> 5 votes
```
# if you have seqtk installed, skip the following two lines
git clone https://github.com/lh3/seqtk
cd seqtk && make
# the main command line
./seqtk cutN -gp10000000 -n1 hg38.fa > hg38-N.bed

```

Option `-n` sets the min stretch length. Just use `-p` as is. It is a bit complicated to explain what it is doing.

# Answer
> 4 votes
This information is stored in the 2bit file representation of the sequence, so if you happen to have a 2bit file locally (or want to download one from UCSC) and have py2bit installed (you'll need version 3.0, since I literally just added support for this):

```
import py2bit
tb = py2bit.open("genome.2bit")
of = open("NNNNN.bed", "w")  # Do change the input and output file names
for chrom in tb.chroms().keys():
    blocks = tb.hardMaskedBlocks(chrom)
    for block in blocks:
        of.write("{}\t{}\t{}\n".format(chrom, block[0], block[1]))
of.close()
tb.close()

```

In case it's helpful, you can also use this to get all soft-masked regions. Just replace `hardMaskedBlocks()` with `softMaskedBlocks()` and ensure to specify `storeMasked=True` when you open the file.

# Answer
> 3 votes
Use twoBitInfo:

```
$ twoBitInfo file.fa -nBed output.bed

```

For example, to get all the N-masked regions on chromosome Y (also note you can use stdout as a filename to write directly to stdout, and use of the url as the input, no need to download the 2bit file):

```
$ twoBitInfo http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.2bit:chrY -nBed stdout | head -3
chrY    0   10000
chrY    44821   94821
chrY    133871  222346

```

You can download twoBitInfo from the directory appropriate to your operating system here: http://hgdownload.soe.ucsc.edu/admin/exe/

# Answer
> 2 votes
I might as well jump on as well. Here's a Perl script that I wrote a while ago to split a fasta sequence at stretches of Ns:

https://github.com/gringer/bioinfscripts/blob/master/fasta-nsplit.pl

I've just modified it to spit out a BED formatted file showing the N locations on standard error. Use it as follows:

```
./fasta-nsplit.pl tig87_withN.fa 2>out.bed > out.split.fa

```

Output (BED file):

```
tig00000087 0   60
tig00000087 900 960
tig00000087 3840    3960
tig00000087 14880   14940
tig00000087 59520   59700
tig00000087 93000   93120
tig00000087 107880  107940
tig00000087 109135  109140

```

Here's the full script, for completeness / fixing / comment:

```
#!/usr/bin/perl
use warnings;
use strict;

my $seq = "";
my $shortSeqID = "";
my $seqID = "";
my $keep = 0;
my $cumLength = 0;
while(<>){
  chomp;
  if(/^>((.+?)( .*?\s*)?)$/){
    my $newID = $1;
    my $newShortID = $2;
    if($seq){
      my $inc = 0;
      while($seq =~ s/(NNNN+)(.*)//){
 my $nStretch = $1;
        my $newSeq = $2;
        printf(">%s.%s\n%s\n", $seqID, $inc++, $seq) if ($seq);
 $cumLength += length($seq);
 printf(STDERR "%s\t%d\t%d\n", $shortSeqID, $cumLength,
        $cumLength + length($nStretch));
 $cumLength += length($nStretch);
        $seq = $newSeq;
      }
      printf(">%s\n%s\n", $seqID, $seq) if ($seq);
    }
    $seq = "";
    $shortSeqID = $newShortID;
    $seqID = $newID;
    $cumLength = 0;
  } else {
    $seq .= $_;
  }
}
if($seq){
  my $inc = 0;
  while($seq =~ s/(NNNN+)(.*)//){
    my $nStretch = $1;
    my $newSeq = $2;
    printf(">%s.%s\n%s\n", $seqID, $inc++, $seq) if ($seq);
    $cumLength += length($seq);
    printf(STDERR "%s\t%d\t%d\n", $shortSeqID, $cumLength,
       $cumLength + length($nStretch));
    $cumLength += length($nStretch);
    $seq = $newSeq;
  }
  printf(">%s\n%s\n", $seqID, $seq) if ($seq);
}

```

# Answer
> 1 votes
Here's a way to generate it yourself from the genome sequence. First, convert the fasta file of the genome to tbl format (`<seq id>\t<sequence>)`, then use perl to find the start and end positions of all stretches of consecutive of `N` or `n`.

```
FastaToTbl hg38.fa |
    perl -F"\t" -ane 'while(/N+/ig){
                            for(0..$#-){
                                print "$F[0]\t$-[$_]\t$+[$_]\n"
                            }
                      }' > hg38.n.bed

```

### Explanation

* `FastaToTbl` : this is a very simple script that converts fasta to tbl. Just save the lines below somewhere in your `$PATH` (e.g. `~/bin`) as `FastaToTbl` and make it executable.

```
#!/usr/bin/awk -f
{
    if (substr($1,1,1)==">")
            if (NR>1)
                printf "\n%s\t", substr($0,2,length($0)-1)
            else
                printf "%s\t", substr($0,2,length($0)-1)
        else
            printf "%s", $0
}
END{printf "\n"}

```
* The Perl magic.

+ The `-a` makes `perl` behave like `awk`, splitting each input line on the value given by `-F` (a tab in this case) and saving the result in the aray `@F`. So, `$F[0]` will be the sequence id.
    + `while(/N+/ig){` : this will match all cases of consecutive `N` or `n` (the `i` flag makes it case-insensitive). Perl will store the start positions of all matches in the special array `@-` and the end positions in `@+`.
    + `for(0..$#-){` : iterate over all numbers from 0 to the final index (`$#-`) in the array `@-`.
    + `print "$F[0]\t$-[$_]\t$+[$_]\n"` : print minimal bed format data. The name of the current sequence (`$F[0]`), the start position of the current match (`$-[$_]`) and the corresponding end position (`$+[$_]`).

I just ran the above on my system and it generated this.

# Answer
> 1 votes
* @user172818 code can't count N at fa begining
* @terdon code get wrong bed position
* i change @terdon code like below; and it works for me

```
awk '{
    if (substr($1,1,1)==">")
        if (NR>1)
            printf "\n%s,", substr($0,2,length($0)-1)
        else
            printf "%s,", substr($0,2,length($0)-1)
    else
        printf "%s", $0
}
END{printf "\n"}' test.fa | \
perl -F"," -ane 'while($F[1] =~ /N+/ig){
                        for(0..$#-){
                            print "$F[0]\t$-[$_]\t$+[$_]\n";
                        }
                  }'

```

---
Tags: genome, bed
---


# Question
Title: Efficiently iterate over bed and subset pyranges objects
I am using pyranges and I have two datasets that I want to compare in several genome-wide intervals based on a bed file.

For every interval in the bed file, I want to get all overlapping positions from the two datasets and sum two columns. That gives me 4 numbers for a fisher exact test, but that's not the most important part.

My current code looks like below, and just iterates over the itertuples() of the bed file as pandas dataframe, and uses that to slice the pyranges object:

```

import pandas as pd
import pyranges as pr
from pyranges import PyRanges
from scipy.stats import fisher_exact
import numpy as np

def main():
    colnames = ["Chromosome", "Start", "End", "num_motifs", "calls", "methylated", "freq", "seq"]
    h1 = PyRanges(pd.read_csv("dataset1.tsv.gz", sep="\t", names=colnames, header=0))
    h2 = PyRanges(pd.read_csv("dataset2.tsv.gz", sep="\t", names=colnames, header=0))
    for chr, begin, end in pr.read_bed("chr21.bed.gz").merge().as_df().itertuples(index=False, name=None):
        odr, pv = fisher_exact([meth_counts(h1[chr, begin:end]),
                                meth_counts(h2[chr, begin:end])])

def meth_counts(interval):
    if not interval.empty:
        all, methylated = interval.as_df()[['calls', "methylated"]].sum()
        return [all - methylated, methylated]
    else:
        return [0, 0]

```

This code works well, but I'm not sure if I use the "full power" of what pyranges could offer me for iterating and subsetting. I uploaded some example data on google drive (250GB).

# Answer
Setup:

```
import pandas as pd
import pyranges as pr
from pyranges import PyRanges
from scipy.stats import fisher_exact
import numpy as np

# ! zcat dataset1.tsv.gz | head -2
# chromosome    start   end num_motifs_in_group called_sites    called_sites_methylated methylated_frequency    group_sequence
# chr21 5010053 5010053 1   3   0   0.000   CACCACGTCCA

# avoid using more columns than we need, especially seq which contains object
# is heavy to lug around
colnames = ["Chromosome", "Start", "End", "calls", "methylated"]
usecols = [0, 1, 2, 4, 5]

h1 = PyRanges(pd.read_csv("dataset1.tsv.gz", sep="\t", names=colnames, header=0, usecols=usecols))
h2 = PyRanges(pd.read_csv("dataset2.tsv.gz", sep="\t", names=colnames, header=0, usecols=usecols))

regions = pr.read_bed("chr21.bed.gz")
regions = regions[~regions.Chromosome.str.contains("_")].merge()

regions.ID = list(range(len(regions)))

regions
# +--------------+-----------+-----------+-----------+
# | Chromosome   | Start     | End       | ID        |
# | (category)   | (int32)   | (int32)   | (int64)   |
# |--------------+-----------+-----------+-----------|
# | chr21        | 5065660   | 5065810   | 0         |
# | chr21        | 5068140   | 5068290   | 1         |
# | chr21        | 5128400   | 5128550   | 2         |
# | chr21        | 5154780   | 5154930   | 3         |
# | ...          | ...       | ...       | ...       |
# | chr21        | 46637040  | 46637190  | 2377      |
# | chr21        | 46638480  | 46638630  | 2378      |
# | chr21        | 46639820  | 46639970  | 2379      |
# | chr21        | 46692680  | 46692830  | 2380      |
# +--------------+-----------+-----------+-----------+
# Unstranded PyRanges object has 2,381 rows and 4 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.

h1
# +--------------+-----------+-----------+-----------+--------------+
# | Chromosome   | Start     | End       | calls     | methylated   |
# | (category)   | (int32)   | (int32)   | (int64)   | (int64)      |
# |--------------+-----------+-----------+-----------+--------------|
# | chr21        | 5010053   | 5010053   | 3         | 0            |
# | chr21        | 5010215   | 5010215   | 5         | 1            |
# | chr21        | 5010331   | 5010335   | 6         | 0            |
# | chr21        | 5010525   | 5010525   | 2         | 2            |
# | ...          | ...       | ...       | ...       | ...          |
# | chr21        | 46699713  | 46699718  | 274       | 254          |
# | chr21        | 46699739  | 46699739  | 133       | 113          |
# | chr21        | 46699756  | 46699756  | 37        | 32           |
# | chr21        | 46699770  | 46699773  | 156       | 130          |
# +--------------+-----------+-----------+-----------+--------------+
# Unstranded PyRanges object has 331,373 rows and 5 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.

```

Get total number of methylated per region:

```
def count_total_and_methylated(df):
    """For each region in bed, get sum of calls and methylated"""

    grpby = df.groupby("ID")
    total = grpby.calls.sum()
    methylated = grpby.methylated.sum()

    chromosome = df.Chromosome.iloc[0]
    start = grpby.Start.first()
    end = grpby.End.first()
    region_id = grpby.ID.first()

    return pd.DataFrame({"Chromosome": chromosome, "Start": start, "End": end, "ID": region_id,
                         "calls": total, "methylated": methylated})

rh1 = regions.join(h1).drop(like="_b").apply(count_total_and_methylated)
rh2 = regions.join(h2).drop(like="_b").apply(count_total_and_methylated)

rh1.File = 1
rh2.File = 2

def add_missing(df, df2):

    """Since join is inner rh1 and rh2 might not have the exact same intervals. Fix this."""

    missing_in_first = np.setdiff1d(df2.ID, df.ID)
    first_missing = df2[df2.ID.isin(missing_in_first)].copy()

    f = df.File.head().iloc[0]

    first_missing.loc[:, "File"] = f
    first_missing.loc[:, "calls"] = 0
    first_missing.loc[:, "methylated"] = 0

    return pd.concat([df, first_missing])

rh1 = rh1.apply_pair(rh2, add_missing)
rh2 = rh2.apply_pair(rh1, add_missing)

j = rh1.join(rh2, suffix=2).drop(like="File*|ID*|(Start|End)2")
j
# +--------------+-----------+-----------+-----------+--------------+-----------+---------------+
# | Chromosome   | Start     | End       | calls     | methylated   | calls2    | methylated2   |
# | (object)     | (int32)   | (int32)   | (int64)   | (int64)      | (int64)   | (int64)       |
# |--------------+-----------+-----------+-----------+--------------+-----------+---------------|
# | chr21        | 5065660   | 5065810   | 188       | 8            | 70        | 0             |
# | chr21        | 5068140   | 5068290   | 168       | 166          | 74        | 72            |
# | chr21        | 5128400   | 5128550   | 901       | 14           | 143       | 1             |
# | chr21        | 5154780   | 5154930   | 653       | 0            | 593       | 0             |
# | ...          | ...       | ...       | ...       | ...          | ...       | ...           |
# | chr21        | 46637040  | 46637190  | 66        | 63           | 69        | 62            |
# | chr21        | 46638480  | 46638630  | 28        | 0            | 30        | 3             |
# | chr21        | 46639820  | 46639970  | 22        | 0            | 18        | 14            |
# | chr21        | 46692680  | 46692830  | 71        | 25           | 57        | 23            |
# +--------------+-----------+-----------+-----------+--------------+-----------+---------------+
# Unstranded PyRanges object has 1,922 rows and 7 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.

```

Compute fisher exact using fisher instead of scipy:

```
def fe(df):
    from fisher import pvalue_npy

    cols = "calls methylated calls2 methylated2".split()
    df[cols] = df[cols].astype(np.uint)
    _, _, two_sided = pvalue_npy(df.calls.values, df.methylated.values, df.calls2.values, df.methylated2.values)

    t1, t2, m1, m2 = df.calls, df.calls2, df.methylated, df.methylated2

    _or =   (t1 * m2) / (t2 * m1)

    df.insert(df.shape[1], "PValue", two_sided)
    df.insert(df.shape[1], "OR", _or)

    return df

result = j.apply(fe, nb_cpu=1)
result.drop(like="methylated*|calls*")
# +--------------+-----------+-----------+------------------------+--------------------+
# | Chromosome   | Start     | End       | PValue                 | OR                 |
# | (object)     | (int32)   | (int32)   | (float64)              | (float64)          |
# |--------------+-----------+-----------+------------------------+--------------------|
# | chr21        | 5065660   | 5065810   | 0.1152954430992567     | 0.0                |
# | chr21        | 5068140   | 5068290   | 0.9999999999995065     | 0.9846955389124064 |
# | chr21        | 5128400   | 5128550   | 0.7075919513120225     | 0.4500499500499501 |
# | chr21        | 5154780   | 5154930   | 1.0                    | nan                |
# | ...          | ...       | ...       | ...                    | ...                |
# | chr21        | 46638480  | 46638630  | 0.24262295081967644    | inf                |
# | chr21        | 46639820  | 46639970  | 0.00024909594620266284 | inf                |
# | chr21        | 46692680  | 46692830  | 0.7354565511951735     | 1.1459649122807019 |
# | chr21        | 36784100  | 36784250  | 1.0                    | nan                |
# +--------------+-----------+-----------+------------------------+--------------------+
# Unstranded PyRanges object has 1,932 rows and 5 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.

```
> 2 votes

# Answer
Unless you need to use Python and can't use `subprocess`, here's a quick CLI one-liner which sums signal from sorted BED5 files, over the genomic space where they overlap:

```
$ bedmap --echo --sum --delim '\t' <(bedops --merge A.bed B.bed ... N.bed) <(bedops --everything A.bed B.bed ... N.bed) > answer.bed

```

(Requires `bash` for process substitutions.)
> 1 votes

---
Tags: python, bed, pyranges
---


title: What are Orthogonal DNA Sequences?

text:
# Question
Title: What are Orthogonal DNA Sequences?
When reading a paper, I came across the following statement

> Standard vectors were constructed, flanked by pairs of homology sequences derived from yeast barcodes at the ends of each expression cassette. We reasoned that as these barcode sequences were designed to be **orthogonal**, they could serve a dual purpose of reducing the probability of mis-annealing and dictating the assembly order of multiple cassettes.

I am not sure what the "orthogonal" means here. I think it is about synthetic biology.

# Answer
> 1 votes
Orthogonal DNA sequences are the shortest possible single stranded nucleotide stretches constructed. These sequences are used as probes in DNA hybridization. Because of its length and construction it has least chance of binding with other probes.

This technique improves the efficacy and specificity and is clinically used to detect genetic disorders in shortest time duration.

Via: https://www.nature.com/articles/srep24650

https://www.pnas.org/content/106/7/2289

---
Tags: phylogenetics
---


title: CDS length for each human gene

text:
# Question
Title: CDS length for each human gene
Does anyone know where and how could I download a list of all human genes and the length of the coding sequence for each gene? Is it possible to do this on the NCBI site, ensembl?

# Answer
While I haven't found a way to limit the results to the canonical transcript only, you can get a list of genes, transcripts and their CDS lengths using Ensemble's BioMart. I have already set it up for you, you can see the results, and modify them, here (click on the "Results" link if you don't see them).

Essentially, you just need to go to BioMart, and

1. select "Ensembl Genes 96" (the number will change if the version changes) as the database and "uman Genes" as the dataset.
2. Click on "Filters", and set `Gene type` to `coding` and `Transcript type` to `protein_coding`.
3. From "Attributes", select whatever you want to see. The "CDS Length" is under "Structures".
> 4 votes

# Answer
Ensembl has an FTP site that allows you to select and download only the coding sequences from many different genomes. https://useast.ensembl.org/info/data/ftp/index.html

To determine the length of those sequences, download the associated gtf or gff3 annotation file. The annotation file is tab delim. The fourth and fifth column represent genomic loci of the annotated region. Subtract the amount in the fourth column from the amount in the fifth column to yield the length of all the annotated features. You can easily do this in R after loading the file in the environment using the Magrittr library. The following code will create a new column called gene\_length with the associated gene lengths.

```
install.packages("magrittr")
# this only needs to be done once
library(magrittr)
# must be run each time the library is neaded
annotation.gtf <- read.table("path/to/annotation.gtf")
annotation.gtf$start <- annotation.gtf[,4]
annotation.gtf$end <- annotation.gtf[,5]
annotation.new-column.gtf <- annotation.gtf %>%
          mutate(gene_length=end-start)

```
> 0 votes

# Answer
Go to the NCBI homepage and search for 'human genome' (the dropdown menu on the left-hand side of the search box should be 'All Databases'). In the result page, click on the 'Download' button; choose 'RefSeq' as the 'Source database' and 'Feature table' as the 'File type' as shown in the screenshot below:

This will download a `genome_assemblies.tar` file which contains the feature table file, `GCF_000001405.39_GRCh38.p13_feature_table.txt.gz`. This file has all of the annotated features in a tab-delimited format with the following fields:

```
              # feature [  1]: gene
                  class [  2]: transcribed_pseudogene
               assembly [  3]: GCF_000001405.39
          assembly_unit [  4]: Primary Assembly
               seq_type [  5]: chromosome
             chromosome [  6]: 1
      genomic_accession [  7]: NC_000001.11
                  start [  8]: 11874
                    end [  9]: 14409
                 strand [ 10]: +
      product_accession [ 11]:
   non-redundant_refseq [ 12]:
      related_accession [ 13]:
                   name [ 14]: DEAD/H-box helicase 11 like 1
                 symbol [ 15]: DDX11L1
                 GeneID [ 16]: 100287102
              locus_tag [ 17]:
feature_interval_length [ 18]: 2536
         product_length [ 19]:
             attributes [ 20]: pseudo

```

You can parse this file to obtain the CDS lengths of all annotated proteins. As an example, you can get a list of all annotated RefSeq proteins and their lengths as follows:

```
zcat GCF_000001405.39_GRCh38.p13_feature_table.txt.gz \
    | awk 'BEGIN{FS="\t";OFS="\t"}($1~/^#/ || $1=="CDS")' \
    | cut -f6-11,13,15,16,18,19

```

The output of the above command will have the following fields:

```
             chromosome [  1]: 1
      genomic_accession [  2]: NC_000001.11
                  start [  3]: 69091
                    end [  4]: 70008
                 strand [  5]: +
      product_accession [  6]: NP_001005484.1
      related_accession [  7]: NM_001005484.1
                 symbol [  8]: OR4F5
                 GeneID [  9]: 79501
feature_interval_length [ 10]: 918
         product_length [ 11]: 305

```
> 0 votes

---
Tags: gene, sequence-analysis, phylogenetics, ensembl
---