# Question
Title: tabix indexing a GTEx variant file- ':' delimited
I have a gtex variant file- the head of which looks as follows:

```
phenotype_id                                    variant_id
chr1:15947:16607:clu_36198:ENSG00000227232.5    chr1_13550_G_A_b38  ...
chr1:15947:16607:clu_36198:ENSG00000227232.5    chr1_14671_G_C_b38  ...
chr1:15947:16607:clu_36198:ENSG00000227232.5    chr1_14677_G_A_b38  ...
chr1:15947:16607:clu_36198:ENSG00000227232.5    chr1_16841_G_T_b38  ...

```

Effectively I would like to be able to lookup the variants inside a particular gene (the file is sorted) and put them in a temporary file: E.g. If the variants in gene "ENSG00000148481- MINDY3" are on line number 87528225 to 87536766- effectively what I would like is the equivalent to `zcat file.gz | sed -n '87528225,87536766p' > MINDY3.txt`. However `zgrep ENSG00000148481 file.gz` is just as fast as the above...

Hence I thought tabix would be the right tool for this-

I would like to tabix index it to make lookups faster. It is gzip compressed and I will firstly have to do:

```
zcat gtex.txt.gz | bgzip > gtex.txt.bgz

```

However I am not quite sure how to proceed from there given that the data is not tab-delimited.

As a trial I tried the first 1000 lines:

```
zcat gtex.txt.gz | head -n 1000 | bgzip > gtex_1000.txt.bgz

./tabix -p bed gtex_1000.gz #index as a bed file
[get_intv] the following line cannot be parsed and skipped: chr1:15947:16607:clu_36198:ENSG00000227232.5        chr1_13550_G_A_b38     ......

./tabix -p vcf gtex_1000.gz #index as a vcf file
Indexing as a bed file results in a warning while indexing as a vcf gives no warning yet either way when I try to retrieve a sequence:

./tabix test.gz chr1:15000:17000

```

It returns nothing.

I am starting to think that I will just have to write a script that splits on the ':' and writes the data to a new file.... and then index that file- which will take a huge amount of time ... Does anyone know of a trick to index the files with unconventional delimiting?

# Answer
> 0 votes
I proceeded as follows:

```
import gzip
import sys

qtls=['sqtls','eqtls']
tissues=['Adipose_Subcutaneous', 'Adipose_Visceral_Omentum', 'Adrenal_Gland', 'Artery_Aorta', 'Artery_Coronary', 'Artery_Tibial', 'Brain_Amygdala', 'Brain_Anterior_cingulate_cortex_BA24', 'Brain_Caudate_basal_ganglia', 'Brain_Cerebellar_Hemisphere', 'Brain_Cerebellum', 'Brain_Cortex', 'Brain_Frontal_Cortex_BA9', 'Brain_Hippocampus', 'Brain_Hypothalamus', 'Brain_Nucleus_accumbens_basal_ganglia', 'Brain_Putamen_basal_ganglia', 'Brain_Spinal_cord_cervical_c-1', 'Brain_Substantia_nigra', 'Breast_Mammary_Tissue', 'Cells_Cultured_fibroblasts', 'Cells_EBV-transformed_lymphocytes', 'Colon_Sigmoid', 'Colon_Transverse', 'Esophagus_Gastroesophageal_Junction', 'Esophagus_Mucosa', 'Esophagus_Muscularis', 'Heart_Atrial_Appendage', 'Heart_Left_Ventricle', 'Kidney_Cortex', 'Liver', 'Lung', 'Minor_Salivary_Gland', 'Muscle_Skeletal', 'Nerve_Tibial', 'Ovary', 'Pancreas', 'Pituitary', 'Prostate', 'Skin_Not_Sun_Exposed_Suprapubic', 'Skin_Sun_Exposed_Lower_leg', 'Small_Intestine_Terminal_Ileum', 'Spleen', 'Stomach', 'Testis', 'Thyroid', 'Uterus', 'Vagina', 'Whole_Blood']

gene_old=''

def index(tissue,qtl,gene_old):
        counter_old,counter=1,1
        if qtl=='sqtls':
                path='/path/to/file/'+tissue+'.v8.sqtl_allpairs.txt.gz'
        else:
                path='/path/to/file/'+tissue+'.allpairs.txt.gz'
        f = gzip.open(path, 'r')
        next(f)
        for line in f:
                try:
                        line=line.decode()
                        line=line[:line.find('\t')]
                        gene=line[line.rfind('ENS'):].split('.')[0]
                        #print(gene,gene_old,counter_old,counter)
                        with open(qtl+'_'+tissue+'.idx','a') as indexed:
                                 indexed.write(str(gene_old)+'\t'+str(counter_old)+'\t'+str(counter)+'\n')
                        gene_old=gene
                        counter +=1
                        counter_old=counter
                except IndexError:
                        print('tissue,qtl,gene_old,counter')
                        counter+=1
        f.close()

for qtl in qtls:
        for tissue in tissues:
               qtl=row.qtls
               tissue=row.tissues
               print(qtl, tissue)
               index(tissue,qtl,gene_old)

```

The above script writes an index for each gene which should look as follows (took 2 days to finish all files):

```
ensid    row_start    row_end
......
ENSG00000230337     1526391 1534329
ENSG00000171819 1534330 1542321
ENSG00000198793    1542322 1550178
ENSG00000120942  1550179 1558010
.....

```

And to get a gene of interest you will have to look up the coordinates for an ensemblid/gene and for example:

```
zcat Brain_Nucleus_accumbens_basal_ganglia.v8.sqtl_allpairs.txt.gz | sed -n '1550179,1558010p;1558011q' > ENSG00000120942.txt

```

# Answer
> 0 votes
Tabix files need a name and a genomic coordinate at least. Tabix can be (ab)used by creating a fake genomic coordinate and then it just indexes names but not sure I would recommend this usage. What is your actual expected application of this?

---
Tags: parsing, tabix, gtex
---


# Question
Title: Compute the significance of the overlap between 2 or more gene sets
I was able to compute the significance of the overlap between 2 gene sets using the cdf function of scipy hypergeometric distribution.

I wish to be able to perform the same calculation for more than 2 gene sets; should I use the multivariate hypergeometric distribution cdf function for that?

Are there any websites that provides the same calculations over gene sets so I can validate my results?

# Answer
This is a shot at it, first an example dataset:

```
import matplotlib.pyplot as plt
import numpy as np
import functools
from matplotlib_venn import venn3

# define universe
uni = ["gene"+str(i) for i in range(1000)]
# some overlap
gs1 = uni[250:300] + uni[900:950]
gs2 = uni[:300]
gs3 = uni[250:500]

```

The ever amazing venn diagram:

```
venn3([set(gs1),set(gs2),set(gs3)],set_labels=["gs1","gs2","gs3"])

```

Then a function to draw a set with length equivalent of each set, randomly from the universe and find length of intersection (all 3):

```
def sim_intersect(uni,set_lengths):
    randomsets = [np.random.choice(uni,n) for n in set_lengths]
    return len(functools.reduce(np.intersect1d,randomsets))

```

We run this 1000 times:

```
permuted_values = [sim_intersect(uni,[len(gs1),len(gs2),len(gs3)]) for i in range(1000)]

plt.hist(permuted_values,bins=range(50))

```

The probability of observing the starting result, using (B+1)/(M+1) as estimator, see this post:

```
(sum(np.array(permuted_values)>obs_n)+1)/(1000+1)
0.000999000999000999

```
> 3 votes

---
Tags: statistics
---


# Question
Title: Pathway level analysis of single-cell gene expression
I'm looking for single-cell specific methods to construct (using gene expression data) new features that express pathway "level" or "activity", and then use these for clustering cells.

One example for bulk RNA-seq is PLAGE, implemented in the GSVA R package.

I have found only one method for single-cell RNA-seq that comes close: pagoda in the scde R package. This performs PCAs for each pathway, using its genes.

A very similar question is Classification (supervised learning) of expression data on pathway level, on bulk/microarray data.

# Answer
> 2 votes
Package Hipathia allows you to compute a value of "activity" for each pathway and each cell, so that the matrix of gene expression is transformed into a matrix of pathway activity. You can use this matrix for further analysis, cell differentiation and clustering. You can also compute functional activity matrices in the same way.

# Answer
> 1 votes
The AUCell R package identifies enriched gene sets within cells. It uses the Area Under the Curve (AUC) to calculate whether a gene set is enriched within the expressed genes for each cell.

# Answer
> 0 votes
Single-Cell\_Signature\_Explorer paper compared several such methods including: AUCell, Seurat CellCycleScore, and GSVA/ssGSEA.

# Answer
> 0 votes
Check our recent paper where we, not only characterized pathway activities in scRNAseq experiments, but also took advantage of mechanistic models in order to make some in silico prediction of drug effects.

---
Tags: r, scrnaseq, pathway
---


title: Core genome too short

text:
# Question
Title: Core genome too short
I'm trying to create a snp distance matrix from `M tuberculosis` isolates to try and infer transmission networks. My plan was to use snippy to make a core genome alignment and then snp-dists to create the matrix but have hit a snag with creating a core genome alignment. When I try it with all 972 of my samples, snippy halts when trying to run snp-sites, saying `Warning: No SNPs were detected so there is nothing to output.` The issue is that with all the samples included, the core genome size is 0.

When looking at the alignment statistics provided by snippy, it looks like some isolates have very low numbers of aligned base pairs (??due to low coverage or contamination). I have progressively filtered out samples with low numbers of aligned base pairs and can then get snippy/snp-sites to work. As I increase my filter stringency, the resulting core genome is still short. Eg Filtering out those with aligned base pairs of \<90% of the reference had a core genome of only 5747bp (from a reference of 4.4 million bp), while excluding 88 (9%) of my samples.

`snippy`'s author recommends using it's output file `core.txt` to figure out which samples are the "bad" outliers. That file provides 1) the length of the reference, 2) the number of aligned base pairs, 3) the number of unaligned base pairs, 3) the number of variant sites, 4) the number of heterogenous sites, 5) the number of masked sites, and 6) the number of low coverage sites.

```
#>   ID                        LENGTH ALIGNED UNALIGNED VARIANT   HET MASKED LOWCOV
#>   <chr>                      <dbl>   <dbl>     <dbl>   <dbl> <dbl>  <dbl>  <dbl>
#> 1 R15795_CATCAAGT_S34_L006 4411532 4192432      6818     755   428 209178   2676
#> 2 R15842_GTCTGTCA_S49_L006 4411532 4187344     10304     784   441 209178   4265
#> 3 R15876_CGCTGATC_S36_L006 4411532 4176662     18992    1292   660 209178   6040
#> 4 R15951_ATTGGCTC_S7_L002  4411532 4170649     14732    1281   980 209178  15993
#> 5 R16019_TGGAACAA_S78_L001 4411532 4190733      8132     712   715 209178   2774
#> 6 R16046_GACTAGTA_S6_L002  4411532 4186069      4309    1257  1084 209178  10892

```

My question is: what heuristic would you use to filter out samples prior to creating a core genome?

And secondarily, what would a reasonable core genome size be for Mtb?

# Answer
Targeting Illumina reads with \>70X depth-of-coverage is a good place to start. Once a robust database is established lower coverage isolates can be used effectively too. More specific to individual SNPs... filtering on VCF QUAL value, Map Quality and/or AC value have shown to be useful.

Core genome will likely be around 1,000 SNPs, but it will depend on the lineage the reference belongs to and the aligned sample lineage. If both are from the same lineage you're likely to be less than 700 SNPs. When from different lineages SNP counts will obviously be higher in addition to the alignment error seen.

Pertaining to filtering individual SNPs...

vSNP is a tool developed for SNP analysis. It can be applied to different specie types but was originally developed for use with TB complex isolates. Code detail is here.

vSNP is a 2-step process. First, generating VCF files. Second, using those VCF files to output SNP tables and trees. The goal of the script is to achieve a high resolution genotyping analysis that can be validated with minimal effort. Some steps are included to help achieve this such as reference selection, automatic grouping and position filtering.

Pseudocode:

For set of VCF files created from the same reference...

```
def find_positions()
    AC=2
    QUAL > 300

```

For each found position...

```
def get_snps()
    AC=2 (solid call)
    QUAL > 50
        THEN --> SNP
    AC=1 (mix call)
    QUAL > 150
        THEN --> AMBIGIOUS
    AC=2
    QUAL >= 50
       THEN --> REF
    QUAL < 50
       THEN --> N
    ALT == "None"
       THEN --> -

```

First find all possible positions that could be informative in all samples/VCF files: `find_positions()`, then define each position to be included in the alignment: `get_snps()`.

Then it can be convenient to filter consistently poor alignment regions such as PPE/PGRS.
> 2 votes

---
Tags: snp
---


# Question
Title: When running ChIPQC bioconductor package after peak calling, should I use sorted BAM files or unsorted BAM files?
I have called peaks using MACS2 and now I'm trying to run ChIPQC bioconductor package. Can someone please tell me should I use sorted BAM files or unsorted BAM files? I have created sorted BAM files and from that indexes for BAM files. When I was trying run ChIPQC earlier, I got an error. Then I figured out that is because I don't have index files. Then I created indexes after sorting BAM files.

# Answer
You need to provide sorted bam files, ideally with an index. If you run `ChIPQCsample`, it calls `ChIPQC:::sampleQC` which iterates through the chromosomes or contigs, using the index, and the absence of the index will crash it.

We can use an example, first in linux:

```
wget ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeHaibTfbs/wgEncodeHaibTfbsA549RxlchV0416101Etoh02AlnRep1.bam -O Etoh02AlnRep1.bam

samtools index ./Etoh02AlnRep1.bam
# s
samtools view -uh Etoh02AlnRep1.bam chr1 > sorted.bam
samtools index sorted.bam
# quick way to unsort
samtools sort -n -o unsorted.bam sorted.bam
# it fails to index
samtools index unsorted.bam
[E::hts_idx_push] Unsorted positions on sequence #1: 246814554 followed by 42922635
samtools index: failed to create index for "unsorted.bam"

```

Then in R:

```
library(ChIPQC)
ChIPQCsample("unsorted.bam")
Bam file has 25 contigs
Error in value[[3L]](cond) :
  failed to open BamFile: failed to load BAM index
  file: unsorted.bam

ChIPQCsample("sorted.bam",annotation=NULL)
Bam file has 25 contigs
Calculating coverage histogram for chr1

Calculating SSD for chr1

Calculating unique positions per strand for chr1

Calculating shift for chr1

 300 / 300
Counting reads in features for chr1

....

```
> 1 votes

---
Tags: r, bioconductor, chip-seq, macs2
---


# Question
Title: DESeq2 - Error in lfcShrink: 'coef' should specify same coefficient as in results 'res'
I have created a DESeq object (dds) and now I want to access the results and apply the lfc shrinkage on them.

```
res <- results(dds, alpha=0.05)
resultsNames(dds)
[1] "Intercept"                                                  "Cell_type_Interstitial_macrophages_vs_Alveolar_macrophages"
[3] "Cell_type_T_cells_vs_Alveolar_macrophages"

```

Here we see all the coefficients I got from the results function. Now I try to use the lfcShrink() function on the second coefficient:

```
res_im_vs_am <- lfcShrink(dds, coef=2, res=res, type="apeglm")
Error in lfcShrink(dds, coef = 2, res = res, type = "apeglm") :
  'coef' should specify same coefficient as in results 'res'

```

it doesn't work, so I try with the name directly copy-pasted from the above code:

```
res_im_vs_am <- lfcShrink(dds, coef = "Cell_type_Interstitial_macrophages_vs_Alveolar_macrophages", res=res, type="apeglm")
Error in lfcShrink(dds, coef = "Cell_type_Interstitial_macrophages_vs_Alveolar_macrophages",  :
  'coef' should specify same coefficient as in results 'res'

```

Same error. According to the `?lfcShrink` for coef: the name or number of the coefficient (LFC) to shrink, consult resultsNames(dds) after running DESeq(dds).

So I do not understand why it is not working. Anyone can help?

# Answer
> 1 votes
We use an example dataset:

```
library(DESeq2)
mat = matrix(rnbinom(9000,mu=100,size=1),ncol=9)
rownames(mat) = paste0("gene",1:1000)
grps = rep(c("A","B","C"),each=3)

dds = DESeqDataSetFromMatrix(mat,data.frame(grps),~grps)
dds = DESeq(dds)

resultsNames(dds)
[1] "Intercept"   "grps_B_vs_A" "grps_C_vs_A"

```

When you call `results`, by default it takes the last coefficient, in this example is C vs A:

```
res <- results(dds, alpha=0.05)

log2 fold change (MLE): grps C vs A
Wald test p-value: grps C vs A
DataFrame with 6 rows and 6 columns
              baseMean     log2FoldChange             lfcSE               stat
             <numeric>          <numeric>         <numeric>          <numeric>
gene1 74.3974631643997 -0.439258650876538  1.22842645044656 -0.357578307367818
gene2 99.4576039995999    1.2903547180366  1.12500005531808   1.14698191518912
gene3 87.1593850381565  -3.11763876685787  1.19848628038255  -2.60131368868297
gene4   51.47821785473  -2.16622171024749 0.953571440023498  -2.27169315200349
gene5 62.4630151830988   1.05879388186169  1.31999477669387  0.802119751195993
gene6 111.967714581788  -1.23916751032406  1.17612862134365  -1.05359863524824

```

If you want to look at the second comparison:

```
res <- results(dds, alpha=0.05,name="grps_B_vs_A")
res_lfc <- lfcShrink(dds, coef=2, res=res, type="apeglm")

log2 fold change (MAP): grps B vs A
Wald test p-value: grps B vs A
DataFrame with 6 rows and 5 columns
              baseMean      log2FoldChange             lfcSE            pvalue
             <numeric>           <numeric>         <numeric>         <numeric>
gene1 74.3974631643997 0.00143892076074213 0.275057978370776 0.987946388824352
gene2 99.4576039995999  -0.031292774512063 0.276439771853046 0.632995722892249
gene3 87.1593850381565 0.00806769808245902 0.274672487627629                NA
gene4   51.47821785473 -0.0394804572502463 0.274244572145665 0.609999235800452
gene5 62.4630151830988   0.065875708916621 0.288411208595865 0.234448669110134
gene6 111.967714581788 -0.0871847461207727 0.295733383995144 0.153081913500716

```

or simply:

```
res_lfc <- lfcShrink(dds, coef=2, type="apeglm")

```

---
Tags: r
---


# Question
Title: How to refine this distribution comparison script
My data look like this

They are compound similarity estimations(the whole file is around 10GB). What I am trying to achieve is to compare the similarity distributions of each compound using the Kolmogorov-Smirnoff test so as to use the data for AIC calculation.

This is my code:

```
import pandas as pd
from scipy.stats import ks_2samp
import re

with open('ecfp4_file.csv', 'r') as f, open('Metrics.tsv', 'a') as f_out:
    f_out.write('compound_1' + '\t' + 'compound_2' + '\t' + 'Similarity' + '\t' + 'KS Distance' + '\n')
    df = pd.read_csv(f, delimiter = ',', lineterminator = '\n', header = 0)

    for m in range(10, 99):
        dc1 = []
        dc2 = []
        m1 = m/100
        m2 = m/100+0.009
        df2 = df.loc[(df['ECFP4'] >= m1) & (df['ECFP4'] < m2)]          #Work on different thresholds each time

        for i in range(0, df2.shape[0]):
            if len(df2.index[df2['id1'] == df2.iloc[i, 0]].tolist()) >= 5:
                dc1.append(df2.iloc[i, 0])
                dc2.append(df2.iloc[i, 1])
        for i in range(0, len(dc1)):
            dc1_l = []
            dc2_l = []
            df3 = df.loc[df['id1'] == dc1[i]]
            dc1_l = df3.iloc[:, 2].tolist()
            df4 = df.loc[df['id1'] == dc2[i]]
            dc2_l = df4.iloc[:, 2].tolist()
            x1 = re.findall(r"statistic=(.*)\,.*$", str(ks_2samp(dc1_l, dc2_l)))
            f_out.write(str(dc1[i]) + '\t' + str(dc2[i]) + '\t' + str(m1) + '\t' + str(x1[0]) + '\n')

```

which so far has taken 10 hours to run and it has not yet passed the 0.1 threshold, so it is highly unpractical. Does anyone know of a better way to do this? Or can help me refine my code for faster and better results?

I also tried to refine it and run it in R, but I get a vector memory exhausted (limit reached) error when I run it

```
install.packages("ggplot2")
library(ggplot2)
PS22_ECFP4 <- read.csv("./PS22_ECFP4.csv")

df <- data.frame(matrix(ncol = 4, nrow = 0))
cols <- c("id1", "id2", "Similarity", "KS Distance")
colnames(df) <- cols

for (i in PS22_ECFP4) {
  sim_df <- PS22_ECFP4[PS22_ECFP4$id1 == i,]
  dist1 <- sim_df[ , 3]
  for (k in 1:nrow(sim_df)) {
    sim <- sim_df[k , 3]
    sim_r <- floor(sim * 100) / 100
    comp2 <- sim_df[k, 2]
    comp2_df <- PS22_ECFP4[PS22_ECFP4$id1==comp2,]
    dist2 <- comp2_df[ , 3]
    if (length(dist1) > 1 & length(dist2) > 1) {
      ks <- ks.test(dist1, dist2)
      df[nrow(df) + 1,] = c(i, comp2, sim_r, ks$statistic)
    }
  }
}

write.csv(df, "./Metrics.csv", row.names = FALSE)

```

# Answer
I am not sure how useful this answer is. But may be you could try this.

For your python script:

You use `f_out.write` within a loop. I think usually this I/O operation is resource intensive. I guess the speed would improve significantly if you store the string in a variable and write the output only once outside the for loop like this:

```
out_str = ''
for i in range(0, len(dc1)):
            dc1_l = []
            dc2_l = []
            df3 = df.loc[df['id1'] == dc1[i]]
            dc1_l = df3.iloc[:, 2].tolist()
            df4 = df.loc[df['id1'] == dc2[i]]
            dc2_l = df4.iloc[:, 2].tolist()
            x1 = re.findall(r"statistic=(.*)\,.*$", str(ks_2samp(dc1_l, dc2_l)))
            out_str +=  (str(dc1[i]) + '\t' + str(dc2[i]) + '\t' + str(m1) + '\t' + str(x1[0]) + '\n')

f_out.write(out_str)

```

Regarding the memory error in your R script, are you using 32bit build of R?

Please see this:

http://web.mit.edu/~r/current/arch/i386\_linux26/lib/R/library/base/html/Memory-limits.html#

They say "Under most 64-bit versions of Windows the limit for a 32-bit build of R is 4Gb: for the oldest ones it is 2Gb. The limit for a 64-bit build of R (imposed by the OS) is 8Tb."

**Conclusion:** Please try using your python script, but use the write operation outside the loop.
> 1 votes

# Answer
If I was doing something like this in R (10G file), I'd use data.table

I can't quite work out what you're R code is doing (does it do what you want it to?)

I think the below might do what you want

```
library(data.table)
dat <- fread("./PS22_ECFP4.csv")
dat[,Similarity:=floor(ECFP4*100)/100]
dat[,'KS Distance':=lapply(1:nrow(dat),function(i)
  ks.test(dat[id1==dat$id1[i],ECFP4],dat[id2==dat$id2[i],ECFP4])$statistic
)]
dat[,ECFP4:=NULL]
fwrite(dat,"./Metrics.csv")

```

Edit: The script loops through the data, and for each row creates two subsets of the data matching id1 and id2. Then does the ks.test on the similarity columns of the two subsets.
> 0 votes

---
Tags: r, python
---


# Question
Title: Transcription factor targets for Human and Mouse genes
I want to get transcription factor targets for my list of genes (For Human and Mouse)

I have found ENCODE Transcription Factor Targets database (The link is given below). But I am looking for some other databases also if available. I will really appreciate if anyone can suggest any other database.

http://amp.pharm.mssm.edu/Harmonizome/dataset/ENCODE+Transcription+Factor+Targets

Thank you

# Answer
JASPAR is another useful transcription factor profiling database.

> The JASPAR CORE contains a curated, non-redundant set of profiles, derived from published and experimentally defined transcription factor binding sites for eukaryotes. It should be used, when seeking models for specific factors or structural classes, or if experimental evidence is paramount.

http://jaspar.genereg.net/
> 2 votes

# Answer
JASPAR is only based on TF motifs, which are very short sequences and not enough to precisely tell you the true TF-motif relationships, especially when multiple TFs have complex logic gates regulating overlapped targets.

An alternative way is to find some transcriptomic data upon perturbation of each TF. The drawbacks are that the redundancy of TFs will mask a lot of targets, and that the up/down regulated genes may not be direct targets but targets of targets (of targets of targets...).

Another way is to find ChIP-seq data for TFs you care about, and hopefully you can find the right biological context. But still ChIP signal isn't exactly binding pattern, which in turn isn't exactly indicating target relationship.

Overall, it's a challenging task to know what are the true targets of a TF.
> 1 votes

---
Tags: genomics, transcriptome, encode
---


# Question
Title: How do you query and explore ENCODE data?
I am looking for a modular way to query data from ENCODE.

For example, I would like to get CHiP-seq or similar tracks for a specific cell line. What's the proper way to do it?

Finally, is there an API to do it?

# Answer
> 7 votes
I don't know whether there is an API, but ENCODE's website does provide an interactive data matrix where you can filter data based on assay and sample type, place data sets in a "shopping cart", and then proceed to "checkout" to download the files of interest.

# Answer
> 5 votes
You can just add `/?format=json` to any page to get the JSON output.

ENCODE REST API documentation: https://www.encodeproject.org/help/rest-api/

Example scripts: https://github.com/ENCODE-DCC/submission\_sample\_scripts

# Answer
> 2 votes
`So let's say I want to get the bed file of this experiment using curl encodeproject.org/experiments/ENCSR000CKC`

You can also search for files that belong to a given experiment, e.g. https://www.encodeproject.org/search/?type=File&dataset=/experiments/ENCSR000CKC/

Then further select file properties (status, format, etc.) using the facets on the left.

This basically builds up a search URL for you, which is also possible to do manually. For example all BED files in ENCSR000CKC with GRCh38 assembly, returned as JSON: https://www.encodeproject.org/search/?type=File&dataset=/experiments/ENCSR000CKC/&file\_format=bed&assembly=GRCh38&format=json (See properties tab here https://www.encodeproject.org/profiles/file.json for all of the file fields you can use.)

You can also specify a `field` parameter to return only fields you are interested in. For example https://www.encodeproject.org/search/?type=File&dataset=/experiments/ENCSR000CKC/&file\_format=bed&assembly=GRCh38&field=s3\_uri&field=cloud\_metadata.url&field=dataset&field=file\_format&format=json only returns the *s3\_uri*, *cloud\_metadata.url*, *dataset*, and *file\_format* fields. These objects are returned in the `@graph` list:

Putting this all together to parse in Python:

```
>>> import requests
>>> import pandas as pd
>>> url = (
...     'https://www.encodeproject.org'
...     '/search/?type=File'
...     '&dataset=/experiments/ENCSR000CKC/'
...     '&file_format=bed'
...     '&assembly=GRCh38'
...     '&field=file_format'
...     '&field=assembly'
...     '&field=s3_uri'
...     '&field=cloud_metadata.url'
...     '&format=json'
... )
>>> r = requests.get(url)
>>> files = r.json()['@graph']
>>> df = pd.DataFrame(files)
>>> df['cloud_metadata'] = df['cloud_metadata'].apply(lambda x: x.get('url'))
>>> df[['@id', 'assembly', 'file_format', 's3_uri', 'cloud_metadata']]
                   @id assembly file_format                                                                                s3_uri                                                                                cloud_metadata
0  /files/ENCFF368QQU/   GRCh38         bed  s3://encode-public/2016/08/04/10493464-917f-4290-bb13-c5fec8319887/ENCFF368QQU.bed.gz  https://encode-public.s3.amazonaws.com/2016/08/04/10493464-917f-4290-bb13-c5fec8319887/ENCFF368QQU.bed.gz
1  /files/ENCFF554NVI/   GRCh38         bed  s3://encode-public/2016/08/04/865f220e-7a47-427a-a7c1-9bdc049d4658/ENCFF554NVI.bed.gz  https://encode-public.s3.amazonaws.com/2016/08/04/865f220e-7a47-427a-a7c1-9bdc049d4658/ENCFF554NVI.bed.gz
2  /files/ENCFF217CRP/   GRCh38         bed  s3://encode-public/2016/08/04/43aa693e-1581-4886-8ace-7e1609f45790/ENCFF217CRP.bed.gz  https://encode-public.s3.amazonaws.com/2016/08/04/43aa693e-1581-4886-8ace-7e1609f45790/ENCFF217CRP.bed.gz

```

The s3\_uri and cloud\_metadata.url both point to the same thing. You can use the s3\_uri to download using awscli or boto3. You can pass the https link to cURL.

# Answer
> 0 votes
There are two places:

> 1. encodeproject.org
> 2. screen.encodeproject.org

Details can be found in the ENCODE3 papers published in July 2020

---
Tags: public-databases, data-download, api, encode, genetics
---


title: Details of ENCODE processing pipelines

text:
# Question
Title: Details of ENCODE processing pipelines
I am trying to find details for how the ENCODE ChIP seq samples were processed.

Linking to a ChIPseq samples on www.encodeproject.org (e.g. https://www.encodeproject.org/experiments/ENCSR051DXE/) offers a PDF "overview" of the pipeline which says basically nothing. The page will also offer you a map of the pipeline, which basically just has headings:

Clicking on a step gives you a popup with a link to the "Chip-seq read mapping" pipeline. This just contains the same overview PDF and a link to "ChIP-seq Data Standards and Processing Pipeline". This page contains again, a vague outline of the pipeline, a link back to the page we just came from and a list of quality control thresholds.

I can't seem to find anywhere the proper details. For example, if BWA was used, was it MEM or aln? What parameters were used for the mapping? One of the steps is read filtering: what where the filters? How were they applied? There is a link to a github page for running the pipeline on DNAnexus, but I can't make head nor tail of it.

# Answer
> 0 votes
The codes are usually available at DCC's GitHub.

You can always email them if you have more questions.

---
Tags: encode, reproducibility
---


# Question
Title: Defining methylation level of CpG island
Differentially methylated site is technically measured by each probe. However, CpG island within gene promoter is comprised of numerous probes. Do we have a consensus of how to determine the methylation level of CpG island using methylation level of probes within corresponding CpG island?

# Answer
> 1 votes
In my experience, the most common approach to find DMR is to use the `DSS` R package to call and merge DML into regions.

DMC/DML (Differential Methylation Cytosine/Loci) is often calculated by comparing the counts of methylated and unmethylated cytosine at one position between samples.

Steps are to first call DMLs with the `callDML()` function. Then call DMRs from these first results with `callDMR()`. Finally annotate the DMLs or DMRs using a `GFF` file. From there you can look for DMR falling into CpG island. Some of the parameters that can be tuned are the **minimum region length**, **minimum number of differential CpG sites**, **percentage of CpG sites required to be significant**. However I have never heard of what would be the standard values for these parameters. It may not be clearly established yet.

So because the process of calling DMLs and DMRs need robust statistical testing, it is suggested to first call the differential methylated states and then annotate it.

---
Tags: methylation
---


# Question
Title: Why doesn't an (Entrez eutils) einfo request for "gene" return the link gene_nucleotide or gene_nucleotide_pos links?
I'm updating a galaxy tool wrapper for Entrez's eutils suite and I'm trying to create a form with valid link selections (among other things) based on the "from" & "to" databases to reduce the size of the select list.

Einfo doesn't return these links for the gene or nuccore databases. I'd found this resource which seemed to show all the possible links, but it it doesn't show these links either. Yet, in working on `elink`'s `acheck` feature, using some gene IDs in a test run, I saw 2 links that were not listed in either resource:

* gene\_nucleotide
* gene\_nucleotide\_pos

I manually checked that I can use those as link names in valid elink queries. It did return valid results (i.e. no error).

So I ran elink on the `gene` database using a few random UIDs to see if I get back a comprehensive list of links (albeit possibly empty results for those links), but it also was missing these 2 links. I'm concerned that I'm not getting all of the possible links for my form interface and that it could produce valid results and that users may want to link via these links, given prior knowledge that they exist.

How do I get a comprehensive list of database links from Entrez and why doesn't einfo return all possible links?

# Answer
I've been corresponding with NLM about this issue and I finally took the time to try out their suggestion (which personally I found hard to see between the lines and which is not a discrete solution, but rather a very time-consuming manual process containing false positives because they say that to get a formal and comprehensive response to my query:

> Collecting the linkname will be a difficult task, that will take time and coordinate/check with relevant groups maintaining individual databases. Your patience will be greatly apprecired. \[sic\]

It took me a bit to decipher their suggestion on how to manually find all links, so I will share what I've learned. There is (currently) no codified or formal means to obtain all possible link names and I infer that many links exist simply because they are utilized in the internal workings of their website. I could be wrong about that, but suffice it to say that there are many undocumented links and they change constantly.

You can get a list of filter items using each individual database's advanced search web interface, "most of which should be the linkname from that source database to other target database". (So this should result in all possible links among a series of false positives.)

Here's how you do that. Let's take the gene database as an example:

1. Go to the advanced search web page for that database, e.g. https://www.ncbi.nlm.nih.gov/gene/advanced
2. In the "Builder" section, select "Filter"
3. Click "Show Index List" to the right of the search term field next to the select list where you selected "Filter" (you will see a multi-select list appear below the search field)
4. Click "Next 200" until you find a selection in that new select list matching the database name followed by "all" and a number in parentheses, e.g. `gene all (29994947)` and click to select that item
5. Click "Refresh index" at the lower right corner next to the multi-select list

The multi-select list will be repopulated and most of the items will be linknames (though you have to replace spaces with underscores in order for those linknames to work with the elink utility).

Doing this for the gene database and then scrolling through just links to nuccore (and its alias "nucleotide" - ignoring est and gss), you will find:

* gene nuccore
* gene nuccore clust
* gene nuccore mgc
* gene nuccore pos
* gene nuccore refseqgene
* gene nuccore refseqrna
* gene nucleotide
* gene nucleotide clust
* gene nucleotide mgc
* gene nucleotide mgc url
* gene nucleotide pos

If you look only at the linknames in the documentation or via einfo (by supplying the `gene` database and looking at similar links to nuccore), you only get these links:

* gene\_nuccore
* gene\_nuccore\_mgc
* gene\_nuccore\_refseqgene
* gene\_nuccore\_refseqrna
* gene\_nucleotide\_mgc\_url

So I believe that given this information and the ever-changing nature of these links, I believe I will allow the user the option to enter a linkname manually, if the set in the select list I generate via einfo does not contain the link they need.
> 3 votes

---
Tags: entrez, e-utilities
---


title: Does NCBI's blast API block my IP?

text:
# Question
Title: Does NCBI's blast API block my IP?
I am trying to run a `blatn` command from both my laptop and within a google colab notebook. I am not sure why, but this command runs properly once and then, on the second execution, it gets stuck for hours without a reply.

Is it possible that NCBI blocks my IP? I've tried also to run `web_blast.pl` and this didn't work either.

biopython:

```
from Bio.Blast import NCBIWWW

fasta = '>seq0\n'+'AACCACCTGACGGGCTCTGT'
blast_res = NCBIWWW.qblast('blastn', 'nt', fasta)
print(blast_res.read())

```

`web_blast.pl` command:

`web_blast.pl blastn nt input.fa > blast.out`

Note, if this is due to throttling by NCBI, I am not sure how do they expect people to query their web service if they block you after one call.

# Answer
You should check the User Guide:

> The NCBI BLAST servers are a shared resource. We give priority to interactive users. In order to ensure availability of the service to the entire community, we may limit searches for some high volume users. Interactive users of the NCBI webpages through a web browser should not encounter problems. We will move searches of users who submit more than 100 searches in a 24 hour period to a slower queue, or, in extreme cases, will block the requests. To avoid problems, API users should comply with the following guidelines:
>
> 1. Do not contact the server more often than once every 10 seconds.
> 2. Do not poll for any single RID more often than once a minute.
> 3. Use the URL parameter email and tool, so that the NCBI can contact you if there is a problem.
> 4. Run scripts weekends or between 9 pm and 5 am Eastern time on weekdays if more than 50 searches will be submitted.
> 1 votes

---
Tags: blast, phylogenetics
---


# Question
Title: How do I write tests for a snakemake pipeline?
I've been playing with snakemake, which is really powerful, however, I'd like to create a 'unit test' for the whole pipeline. e.g. a way to call the pipeline that tests that all the rules work together, all the steps complete, etc.

I assumed I could find some guides about how to do this, but I don't find any clear information about how to tests pipelines.

Any advice or suggestions would be welcome.

# Answer
> 5 votes
You can run CI on a snakemake workflow when provoding a minimal example dataset. Maybe an example from the snakemake-workflows repo will be helpful: https://github.com/snakemake-workflows/dna-seq-varlociraptor

# Answer
> 2 votes
This is what I do when I want to check whether a program runs from start to finish and produces the expected output on test data. Here I used it to test a C program but the idea can be applied also to a snakemake pipeline as in your case.

The idea is pretty simple really: use python's `subprocess` module to run the pipeline then check the results using the `unittest` module (or your favourite unit test package).

For example, save this in `test.py`:

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

import unittest
import os
import subprocess as sp

class Tester(unittest.TestCase):

    def testPipelineMakesStuff(self):
        # Run the pipeline
        p= sp.Popen('snakemake [some parameters]', shell=True, stdout= sp.PIPE, stderr= sp.PIPE)

        # Optionally, get stdout and stderr
        stdout, stderr= p.communicate()

        # Check exits code and other expected output
        self.assertEqual(0, p.returncode)
        self.assertTrue(os.path.isfile('some.pipeline.output.txt'))

        ...Maybe read some.pipeline.output.txt and check it looks ok

    def testPipelineReturnsWithError(self):
        p= sp.Popen('snakemake [some parameters]', shell=True, stdout= sp.PIPE, stderr= sp.PIPE)
        stdout, stderr= p.communicate()
        self.assertTrue(p.returncode != 0)
        self.assertTrue('Something went wrong' in stderr.decode())

    def ...more test cases...

if __name__ == '__main__':
    unittest.main()

```

Run this script as:

```
 python test.py

```

and check all tests pass ok. Of course, you can embed `python test.py` in travis or other continuous integration system as I've done here.

---
Tags: snakemake, reproducibility
---


title: Inheritance annotation for nsSNPs

text:
# Question
Title: Inheritance annotation for nsSNPs
I am relatively new to SNP analysis. Is there a database to find annotation about the Mendelian inheritance of SNPs?

I have a small list of nsSNPs (non-synonymous SNPs) and I need to find how they are they passed to from parents to offspring.

Thank you in advance.

# Answer
> 2 votes
If you only have a small number of genes, OMIM is a very easy way to do this.

For example:

If you have too many to do manually, I suspect there is an easy way to query lists, although I'm not sure exactly the way myself. Perhaps someone can add a new answer if so.

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


# Question
Title: Error using htseq-count: Could not retrieve index file
I have a total RNAseq dataset that I aligned using STAR producing BAM files (sorted by coordinates). I am now trying to get counts for the lncRNA sequences using htseq-count, with the command:

```
htseq-count -m intersection-nonempty -f bam -i gene_id -r pos -s reverse ./OVCA429_Aligned.sortedByCoord.out.bam ./GRCh38_latest_genomic.gtf > ./OVCA429.counts.txt

```

This is the output:

```
[E::idx_find_and_load] Could not retrieve index file for './OVCA429_587680_1_1.rnaseq.fastqAligned.sortedByCoord.out.bam'
100000 GFF lines processed.
200000 GFF lines processed.
300000 GFF lines processed.
400000 GFF lines processed.
...(truncated output here)...
3775001 GFF lines processed.
[E::idx_find_and_load] Could not retrieve index file for './OVCA429_587680_1_1.rnaseq.fastqAligned.sortedByCoord.out.bam'
100000 alignment records processed.
200000 alignment records processed.
300000 alignment records processed.
...(truncated output here)...

Output text file:
__no_feature    20368875
__ambiguous 0
__too_low_aQual 0
__not_aligned   0
__alignment_not_unique  2634439

```

I've read through the htseq-count literature but it never asks to include an index file, therefore i'm not sure how to deal with this problem. Any help/ suggestions would be welcome! Thank you, in advance.

# Answer
> 0 votes
The problem seems to be missing index files of the bam files (therefore the error message that the index file is missing). You should be able to generate those with `samtools`, e.g.

```
samtools index ./OVCA429_587680_1_1.rnaseq.fastqAligned.sortedByCoord.out.bam

```

---
Tags: rna-seq, sam, errors, featurecounts
---


# Question
Title: How to find possible viral subsequences in a Eukaryote genome?
I have an assembly of an algae and want to find if it has anything analogous to endogenous retroviruses' subsequences in its DNA.

My guess is I should manually blastn all subsequences until I find (or not something) that could have come from a virus.

Is there a better way to do this?

# Answer
Your guess sounds reasonable to me but you should consider using KRAKEN2. You can create your custom database (you would have done it with makeblastdb anyway) and kraken will check k-mers against the database to find homologous sequences.
> 2 votes

# Answer
Searching directly for a reverse transcriptase may find recently incorporated retroviruses.
> -1 votes

---
Tags: blast, blastn
---


title: RCSB API: limited returned result

text:
# Question
Title: RCSB API: limited returned result
Is there a limit of query results on RCSB API?

Interestingly, no matter how to adjust my query criteria, only 10 assemblies are returned. Here is an example query json: (truncated from a longer query file)

```
{
  "query": {
    "type": "group",
    "logical_operator": "and",
    "nodes": [
      {
        "type": "terminal",
        "service": "text",
        "parameters": {
          "operator": "greater",
          "value": 1,
          "attribute": "rcsb_assembly_info.polymer_entity_instance_count"
        }
      }
    ]
  },
  "return_type": "assembly"
}

```

Example python code requesting query:

```
url_query = 'https://search.rcsb.org/rcsbsearch/v1/query?json='
query_str = open(query_json, 'r').read()
query_str = query_str.replace('\n', '').replace(' ', '')
url = url_query + query_str

r = requests.get(url)
content = r.content.decode('utf-8')
result = json.loads(content)
print(len(result['result_set']) # always equal to 10

```

# Answer
> 2 votes
Staff from rcsb kindly reminds me the solution \- add the following lines in the json query.

```
"request_options": {
"return_all_hits": true
}

```

---
Tags: assembly, api
---


# Question
Title: Odds ratio and enrichment of SNPs in gene regions?
I did a QTL analysis with a panel of 7M SNPs, and want to analyze the enrichment of the significant qtl-SNPs in different genic regions (promoters, gene bodies, TFBS, etc.).

A straightforward way to do it would be with an odds ratio (OR), using as 'background' the full 7M SNPs, but I prefer to compare with a background with a similar MAF distribution to that of the qtl-SNPs. For that, I'm thinking in two alternatives:

1. Take a single **random sample of all SNPs** with size 10*n* (*n* is the number of qtl-SNPs), keeping the same MAF distribution of the qtl-SNPs results, as well as the same proportion of qtl-SNPs in the sample compared to the full 7M panel. Then, test enrichments using ONLY the random sample (with OR and Fisher's exact test).
2. Do a **bootstrapping**, taking multiple random samples (excluding qtl-SNPs), each of size *n*, and keeping as well the same MAF distribution. For each iteration, get the OR between the full qtl-SNPs and the random sample. Then, the mean of the OR will be the point estimate and the quantiles the confidence intervals. For this approach, the problem would be to compute a *p*-value.

Which approach seems to fit better? If the second one, how to quantify the significance? Can you suggest a better approach?

Thanks.

# Answer
> 3 votes
The first approach only address the question of how likely are you to end up with the observed over-representation given the MAF distribution.

My suggestion is to use the second approach, but I am not sure if you would call it bootstrap. Bootstrap in general means sampling with replacement to estimate the uncertainty of a parameter, in this case, OR. So even if you resample your full dataset, you get an estimate of how variable your enrichment is, not how likely you are to observe it by chance. so just name it permutation or sampling.

This is what I think might work

1. Calculate your true OR, qtl snps versus all others for say TFBS
2. Divide your snps according to MAF into 10 bins / stratas
3. Within each bin, we permute the label whether the snp is a qtl
4. Calculate OR of selection / others for TFBS
5. Repeat 1-3 N times to collect the permutated OR
6. p value would be number of (permutated OR \> observed OR+1)/(N+1)

Below is an example implementation in R, first we have 1e5 signficant snps and they are 5 times as likely to be in TFBS:

```
snp_is_qtl = rep(c(1,0),c(0.1e6,6.9e6))
snp_is_TFBS = c(sample(0:1,0.1e6,replace=TRUE,prob=c(0.95,0.05)),
sample(0:1,6.9e6,replace=TRUE,prob=c(0.99,0.01)))

```

We have the MAF for all snps, with those in qtl being different from others and put them in bins.

```
snp_MAF = c(runif(0.1e6,0.4,0.8),runif(6.9e6,0.1,0.9))

```

We sort the snps according to MAF to make it easier:

```
o = order(snp_MAF)
snp_is_qtl = snp_is_qtl[o]
snp_is_TFBS = snp_is_TFBS[o]
snp_MAF = snp_MAF[o]

```

And we create the bin:

```
bin_MAF = cut(snp_MAF,10)

```

Lastly a function to calculate OR:

```
calOR = function(snps,TBFS){
tab = table(snps,TBFS)
log(tab[1,1])+log(tab[2,2]) -log(tab[1,2]) -log(tab[2,1])
}

```

So our observed OR (in log) is:

```
obs_OR = calOR(snp_is_qtl,snp_is_TFBS)
obs_OR
[1] 1.656859

```

Now we collect the OR under permutation:

```
N = 100
permutated_OR = sapply(1:100,function(i){
print(i)
perm_qtl = unlist(tapply(snp_is_qtl,bin_MAF,sample))
calOR(perm_qtl,snp_is_TFBS)
})

hist(permutated_OR)

```

Your observed OR way exceeds the permuted, even after controlling for MAF. Note the computation is a bit heavy because of the large number of snps. So it makes sense to calculate all the ORs (i.e TFBS, promoters..) with each permutation

---
Tags: statistics, snp, gwas, p-values, qtls
---


title: Setting Contrast for DESeq2 results

text:
# Question
Title: Setting Contrast for DESeq2 results
*This question was also asked on BioStars*

I am trying to recreate this heatmap. How would I compare the four variables: Ly49+/- and MOG/MOGSP to get a single heatmap? Thank you for helping me through this.

I need to compare both Ly49+ versus Ly49− and MOG versus MOG plus SP. Column data in at the bottom. The function I used for DESeq2 is,

```
dds = DESeqDataSetFromMatrix(countData=countData,

colData=colData,

design= ~ cell_treatment)

```

I am not sure if the design is right. For the results, how should I write the contrasts? I am very confused and would really appreciate any help!

```
DataFrame with 11 rows and 2 columns

cell treatment

<character> <character>

MOGSP_2 SPtetramer+CD8+TCell MOGSP

MOGSP_3 SPtetramer+CD8+TCell MOGSP

MOG_3 SPtetramer+CD8+TCell MOG

MOG_4 SPtetramer+CD8+TCell MOG

MOG_5 SPtetramer+CD8+TCell MOG

Ly49N_1 Ly49-CD8+TCell MOGSP

Ly49P_1 Ly49+CD8+TCell MOGSP

Ly49N_2 Ly49-CD8+TCell MOGSP

Ly49P_2 Ly49+CD8+TCell MOGSP

Ly49N_3 Ly49-CD8+TCell MOGSP

Ly49P_3 Ly49+CD8+TCell MOGSP

```

Thank you!!

# Answer
Your table is not very clear. If I understood it correctly, there are three groups of t-cells, SPtetramer+, Ly49- and Ly49+ and within it there is a MOGSP vs MOG comparison. I would set colData like this:

```
                   name             Ly treatment
1  SPtetramer+CD8+TCell SPtetramerplus     MOGSP
2  SPtetramer+CD8+TCell SPtetramerplus     MOGSP
3  SPtetramer+CD8+TCell SPtetramerplus       MOG
4  SPtetramer+CD8+TCell SPtetramerplus       MOG
5  SPtetramer+CD8+TCell SPtetramerplus       MOG
6        Ly49-CD8+TCell      Ly49minus     MOGSP
7        Ly49+CD8+TCell       Ly49plus     MOGSP
8        Ly49-CD8+TCell      Ly49minus     MOGSP
9        Ly49+CD8+TCell       Ly49plus     MOGSP
10       Ly49-CD8+TCell      Ly49minus     MOGSP
11       Ly49+CD8+TCell       Ly49plus     MOGSP

```

I generate an example dataset:

```
counts = matrix(rnbinom(11000,mu=100,size=1),ncol=11)

```

And in the design you specify these two factors, meaning gene expression can be affected by these two conditions:

```
dds = DESeqDataSetFromMatrix(counts,colData=colData,~Ly +treatment)
dds = DESeq(dds)

```

You can compare them like this:

```
head(results(dds,c("Ly","Ly49plus","Ly49minus")))
log2 fold change (MLE): Ly Ly49plus vs Ly49minus
Wald test p-value: Ly Ly49plus vs Ly49minus
DataFrame with 6 rows and 6 columns
          baseMean     log2FoldChange            lfcSE               stat
         <numeric>          <numeric>        <numeric>          <numeric>
1 91.4458114779915  0.712987165797688 1.28291188260395  0.555756927241586
2 81.7945234843954 -0.418266530066077 1.14378290088875 -0.365686993345568
3 94.4351681904941 -0.775592496878288 1.32965355129901 -0.583304196886903
4 107.484335203298  -1.36694885973712 1.38595911743023 -0.986283680770933
5 61.2205696982841 -0.760219864564472 1.16779888688769 -0.650985262188887
6 95.5908168437137  0.497427543096335 1.12768919469662   0.44110340458672
             pvalue              padj
          <numeric>         <numeric>
1 0.578377034355333 0.977537706056193
2    0.714598652706 0.992321035598301
3 0.559688537757548 0.977537706056193
4 0.323993925580938 0.954772990635094
5 0.515056000552832 0.977537706056193
6 0.659138138867294 0.992321035598301

```

To compare groups within another factor:

```
results(dds,c("treatment","MOGSP","MOG"))
log2 fold change (MLE): treatment MOGSP vs MOG
Wald test p-value: treatment MOGSP vs MOG
DataFrame with 6 rows and 6 columns
          baseMean     log2FoldChange            lfcSE               stat
         <numeric>          <numeric>        <numeric>          <numeric>
1 91.4458114779915 -0.875905107403943 1.43007567177694 -0.612488642867122
2 81.7945234843954  0.595588650280833 1.28109080878515  0.464907441530728
3 94.4351681904941 0.0255578400583986 1.48592401497315 0.0171999643325371
4 107.484335203298  -1.46652057431159 1.55285829334496 -0.944400774105796
5 61.2205696982841  -1.27442596494316 1.30489822606005  -0.97664778715433
6 95.5908168437137  -2.35402796630456 1.26608624796739  -1.85929510733078
              pvalue              padj
           <numeric>         <numeric>
1  0.540214509839105 0.944881849859735
2  0.641997741235223 0.973651054534338
3  0.986277090644156 0.994944323499595
4  0.344964886056295  0.87421063696937
5  0.328743552427222                NA
6 0.0629853201285643 0.525889480109579

```

The colData i used:

```
structure(list(name = structure(c(3L, 3L, 3L, 3L, 3L, 1L, 2L,
1L, 2L, 1L, 2L), .Label = c("Ly49-CD8+TCell", "Ly49+CD8+TCell",
"SPtetramer+CD8+TCell"), class = "factor"), Ly = c("SPtetramerplus",
"SPtetramerplus", "SPtetramerplus", "SPtetramerplus", "SPtetramerplus",
"Ly49minus", "Ly49plus", "Ly49minus", "Ly49plus", "Ly49minus",
"Ly49plus"), treatment = structure(c(2L, 2L, 1L, 1L, 1L, 2L,
2L, 2L, 2L, 2L, 2L), .Label = c("MOG", "MOGSP"), class = "factor")), class = "data.frame", row.names = c(NA,
-11L))

```
> 3 votes

---
Tags: rna-seq, deseq2, design
---


# Question
Title: Pointing Cromwell to local EBS mounted storage on AWS Batch
I'm trying to run the WDL tasks from the GTEx RNA-Seq pipeline with Cromwell using AWS Batch as a backend.

I store the STAR alignment index on an Elastic Block Store since my Cromwell server instance lacks sufficient disk space. The EBS Block is mounted to my Cromwell server at `/data`. The STAR index is located at `/data/genomic-references/rnaseq/star_index_oh50`.

Try as I might, I cannot get Cromwell to run the STAR task correctly.

My JSON inputs file:

```
 {
   "rnaseq_pipeline_fastq_workflow.prefix": "s3://mybucket/test/CL100148068",
   "rnaseq_pipeline_fastq_workflow.star.fastq1": "s3://anotherbucket/genotypedata/sample1.fq.gz",
   "rnaseq_pipeline_fastq_workflow.star.star_index": "/data/genomic-references/rnaseq/star_index_oh50/",
   "rnaseq_pipeline_fastq_workflow.star.memory": "32",
   "rnaseq_pipeline_fastq_workflow.star.num_threads": "4",
   "rnaseq_pipeline_fastq_workflow.star.disk_space": "200",
   "rnaseq_pipeline_fastq_workflow.star.num_preempt": "0"
 }

```

Error message:

```
tar: /data/genomic-references/rnaseq/star_index_oh50/: Cannot open: No such file or directory

```

Note that `/data/genomic-references/rnaseq/star_index_oh50/` does exist and is readable.

The relevant command in the GTEx container is here:

```
tar -xvvf ${star_index} -C star_index --strip-components=1

```

If the files are stored locally, then why does `tar` fail to wrap the directory?

# Answer
Somebody helped me with this one through other channels, but I wanted to make it publicly available:

Make a tarball of the STAR index. The option `--strip-components=1` strips the first level of the path from the untarred items; thus, the expected tarball was made with something like

```
tar -cvf ${star_index} the_star_index_directory/

```

Then point Cromwell/WDL to the tarball, e.g. on S3. The GTEx container will unwrap the tarball on the fly during execution.
> 0 votes

---
Tags: rna-seq, genomics, wdl, gtex, aws
---


# Question
Title: Details of DESeq2 modeling a batch effect
When correcting my data for a batch effect using removeBatchEffect, some of the gene expression values become negative.

When searching for differentially expressed genes, I do not use the data above, but rather model the batch using deseq2 (design=~Batch + Condition).

However, I started worrying. When DESeq2 introduces a batch into the model, does it allow for negative values "behind the scenes"?

If it does, I do not understand how that can make sense, in the context of RNA-Seq data.

# Answer
> 1 votes
DESeq2 uses the batch information (and everything else in the design) to produce offsets for its GLM. For a background on that please check how linear models work, e.g. using the StatQuest series of statistics videos over at YouTube. It still operates on the raw counts. The same goes for the normalization factors. `removeBatchEffect` fits a linear model to the data including the batch information and then subtracts the batch component from the counts (that is basically the baseline difference). If you are interested in preserving the integer nature of counts and preserving zeros as smallest values after explicit batch correction you may want to check `ComBat-seq()` from the sva package. It operates on raw counts and returns batch-corrected raw counts which you could then normalize and calculate CPMs from (if you need batch corrected CPMs with no negative values). I find it useful and prefer it over `removeBatchEffect` as it avoids the unfortunate negative counts which sometimes messes up plotting scripts that expect zero as the smallest possible value.

---
Tags: deseq2, limma, batch-effects
---


# Question
Title: How is the GT field in a VCF file defined?
As my question in SO was closed and asked to be posted in this forum, I am posting it here.

I am not from the bioinformatics domain. However, for the sake of analysis, I am trying to pick up certain basics related to the GT field in the VCF file.

I know we have a field called `Alleles`. May I know under what circumstances GT takes a certain value and how they are called? Can you confirm my understanding?

```
Ref Alt    GT   Name
A    A     0/0  Homozygous
A    G     0/1  Heterozygous (does 1/0 also mean the same?) What's this 0 and 1 actually?
A   [C,CA] ??    ??
??   ??    1/1  HOM_ALT? How and why?

```

Can the experts here help me to fill the `question marks` and also help me understand with the Ref and Alt combinations when a genotype can take a value of 0/0 or 1/0 or 1/1 or 1/2 etc. and what are the names for those values? Like when it is called HOM\_ALT etc.?

Any simple explanation for beginner like me (with no background in bioinformatics/biology) can be helpful.

# Answer
You can get most of the info from this paper. See Fig. 1 and the surrounding text. Quoting from there, "GT, genotype, encodes alleles as numbers: 0 for the reference allele, 1 for the first allele listed in ALT column, 2 for the second allele listed in ALT and so on."

In your case, the reference allele, here a single nucleotide, is A. When the alternate allele is also A, the genotype GT is reference, or 0. There are 2 copies of each allele in the human genome in non-sex chromosomes chr1-chr22, hence 0/0, or homozygous reference or HOM\_REF.

When ALT=G, and the GT column is 0/1, this means that you have 1 reference allele (0), and 1 alternate allele (1). This means that you have A on one copy of this locus, and G on another. The convention is write GT field in ascending order, so 0/1 rather than 1/0. This is called heterozygous, or HET.

When ALT=C,CA, the GT is probably 1/2, because there are 2 alternate alleles, and I assume we continue with the same chromosome present in 2 copies. This means there are no reference alleles here at all, only alternate alleles. It is a heterozygous genotype composed of two different ALT alleles, or HET\_ALT. Note that it is not enclosed in square brackets in the vcf file format: `A <tab> C,CA <tab> 1/2 ...`.

Finally, these are some examples of HOM\_ALT:

```
A    C     1/1
A    G     1/1
A    CA    1/1

```

This means that the same ALT allele (either C, or G, or CA) is present in 2 copies. There is no reference allele present. This is called homozygous alternate genotype.

In general, the name *homo* means the same, and *hetero* means different, in the context of genotypes.

**REFERENCES:**

Danecek P, Auton A, Abecasis G, et al. The variant call format and VCFtools. Bioinformatics. 2011;27(15):2156-2158. doi:10.1093/bioinformatics/btr330 : https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3137218/

**SEE ALSO:**
> 7 votes

# Answer
0 and 1 is just a way of coding the reference (0) and alternate (1) allele. These could be A/G, C/T etc. It's just a simplified way of expressing the different alleles.

`0/1` and `1/0` functionally mean the same thing (that the individual is a heterozygote) - since the genotype is unphased, the alleles aren't ordered. The `/` symbol tells you the genotype is unphased. However, in practice, a heterozygous genotype is always written as `0/1` as a matter of convention. `0/0` is referred to as homozygous reference, and `1/1` as homozygous alternative.

`0|1` and `1|0` **do** mean something different, since the pipe symbol `|` tell us the order of the alleles matters.

Occasionally you will have multi allelic positions, where you have more than one alternate allele, and therefore the field will look like:

```
Ref   Alt
A   [C,GT]

```

The extra alt allele is given the number 2, so in this case if an individual had the genotype CT then their genotype code would be 0/2.

Full detail is given in the official VCF specification page.
> 7 votes

---
Tags: vcf, file-formats
---


title: BioJava - Protein to DNA

text:
# Question
Title: BioJava - Protein to DNA
Using Biojava 5, I am trying to convert a sequence of Amino Acids to a Nucleotide sequence. Any idea how to do this?

# Answer
> 1 votes
Note that unlike the reverse (nucleotide sequence to protein sequence), which is uniquely defined given a reading frame, the mapping from protein sequence to nucleotide sequence is not unique (it is not a function in the mathematical sense). The method you use to select which of the many (exponential in input sequence length) possible sequences to output depends critically on what the question is.

If the goal is to design a synthetic gene that will produce a desired protein, then you need to optimize it for expression in the organism into which the gene will be inserted, subject to the constraint of the encoded amino acid sequence.

If the goal is to search for nucleotide sequences in one or more genomes that could give rise to a certain polypeptide, there is a much better way than to do reverse translation at all. In particular, most putative protein-coding open reading frames in most sequenced genomes have already been run through forward translation code and the expected amino acid sequence encoded in them can be found in databases. So, this problem is reduced to a substring search in these databases.

---
Tags: phylogenetics
---


title: Get Dockerfile of a bioconda package

text:
# Question
Title: Get Dockerfile of a bioconda package
Bioconda packages have associated containers images. E.g. `bioconductor-csaw` has the image `quay.io/biocontainers/bioconductor-csaw`.

I would like to customize the image, in particular, I would like to add rstudio.

It is possible to get the Dockerfile associated to the image of a bioconda package?

How do you suggest to proceed to customize the image?

# Answer
> 1 votes
The Dockerfile is generated by galaxy-utils, so I'm not entirely sure what it looks like. However you should be able to generate such a "mulled" container as documented here. In short, something of the form:

```
mulled-build build 'bioconductor-csaw=someversion--somebuildnumber' 'someotherpackage=someversion--somebuildnumber' ...

```

That probably requires setting the CONDA\_IMAGE environment variable. Bioconda currently sets that to `quay.io/dpryan79/mulled_container:latest`.

---
Tags: rstudio, bioconda, docker
---


title: TopHat2 versus HISAT2 inner workings

text:
# Question
Title: TopHat2 versus HISAT2 inner workings
In my intro to bioinformatics course, we mentioned that TopHat2 and HISAT2 will both try to align as many reads as possible to the reference genome (TopHat2 has been superseded by HISAT2). For the reads that were not mapped, these probably cross exon boundaries. To solve this issue, both TopHat2 and HISAT2 chop up the unmapped reads into around 25bp fragments. Then: "TopHat2 tries to align many fragments from each read using all possible exons as the search space. HISAT2 focuses on one fragment from each read, and once that fragment has been anchored it only searches locally."

My question is about that last part "once that fragment has been anchored, \[HISAT2\] only searches locally." Isn't it possible that a given fragment will map to many different regions of the genome? So how is it possible that it will be "anchored"? Are we going to anchor it multiple times and see which anchor makes most sense or what exactly?

# Answer
I don't believe this statement is true:

> To solve this issue, both TopHat2 and HISAT2 chop up the unmapped reads into around 25bp fragments.

TopHat2 did not support spliced alignment directly. This is why it followed this general approach of chopping the read up into small parts, followed by "contiguous" alignment of the sub-reads. Hisat2 / Hisat-genotype support spliced alignment directly insofar as they don't pre-partition the reads into small parts, but allow dynamic extensions of alignments across splice junctions (in a way similar to how STAR works). This is generally a fundamentally better approach and is one reason, among many, why TopHat/TopHat2 have been put into maintenance-only mode and have been officially superseded by Hisat2/genotype.
> 2 votes

---
Tags: sequence-alignment
---


title: get sequences begining with TA

text:
# Question
Title: get sequences begining with TA
I have a fasta file and the sequences in them are arranged like this:

```
>chr8:106028475-106028517()
GTaccaggcctagggttgcactcctgtagtctcagcattggg
>chr8:106040430-106040472()
TAACCGGCTCTCCGCTGCCCTCATCTGTCTAGAATTCAGGTA
>chr9:100245605-100245648()
TACCACCTCCCTCACTCTTTCCCAAATTTATGAGTTAGGCTAT

```

I want to extract only those which begin with TA. so my final output should look like

```
>chr8:106040430-106040472()
TAACCGGCTCTCCGCTGCCCTCATCTGTCTAGAATTCAGGTA
>chr9:100245605-100245648()
TACCACCTCCCTCACTCTTTCCCAAATTTATGAGTTAGGCTAT

```

Kindly help.

# Answer
> 1 votes
seqmagick has a handy `--seq-pattern-include` feature that works similarly to what Kamil did with a regular grep.

```
$ cat sequences.fasta
>chr8:106028475-106028517()
GTaccaggcctagggttgcactcctgtagtctcagcattggg
>chr8:106040430-106040472()
TAACCGGCTCTCCGCTGCCCTCATCTGTCTAGAATTCAGGTA
>chr9:100245605-100245648()
TACCACCTCCCTCACTCTTTCCCAAATTTATGAGTTAGGCTAT

$ seqmagick convert --seq-pattern-include "^TA" sequences.fasta -
>chr8:106040430-106040472()
TAACCGGCTCTCCGCTGCCCTCATCTGTCTAGAATTCAGGTA
>chr9:100245605-100245648()
TACCACCTCCCTCACTCTTTCCCAAATTTATGAGTTAGGCTAT

```

It's a great little tool for this sort of thing built around Biopython, though it's been falling behind in compatibility with the latest versions as you can see from some of the open GitHub issues.

# Answer
> 0 votes
Assuming your sequences are always on a single line

```
grep -B 1 "^TA" seqences.fasta | grep -v "\-\-" > filtered_seqences.fasta

```

should do the job. Argument `-B` specifies how many lines before the pattern match should be printed. But there grep separates individual entries with `--` that's why I used a second grep to get rid of those.

Ha, apparently there is an undocumented `--no-group-separator` flag, so this should do as well

```
grep -B 1 "^TA" seqences.fasta --no-group-separator > filtered_seqences.fasta

```

---
Tags: fasta, sequence
---


# Question
Title: Extracting all information about a sample when using xtract from e-utilities
I would like to extract all information about each `SAMPLE` after running the following query (run the query and add a `| grep SAMPLE` for clarification on what I mean by `SAMPLE`):

`esearch -db sra -query PRJNA514750 | efetch -format xml`

I tried the following:

`esearch -db sra -query PRJNA514750 | efetch -format xml | xtract -pattern EXPERIMENT -element SAMPLE`

but this returns nothing (PS: `SAMPLE`s are within an `EXPERIMENT` tag). I read in the e-utilities guide that `-pattern` will divide the data into rows, and `-element` into columns, so I'm presuming that this didn't work because `SAMPLE` has multiple tags within it. So I then tried:

```
esearch -db sra -query PRJNA514750 | efetch -format xml | xtract -pattern SAMPLE -element random_SAMPLE_tag

```

where `random_SAMPLE_tag` is any tag within `SAMPLE`.

Here's a concrete example:

```
esearch -db sra -query PRJNA514750 | efetch -format xml | xtract -pattern SAMPLE -element TITLE

```

This works, but I want to get all the information about each `SAMPLE`, and I do not know beforehand what the tags within it are (I manually got `TITLE` in this case), and since I want to get this info for a quite a few studies, I can't manually check this.

# Answer
You can give pysradb a shot.

```
study_accession experiment_accession experiment_title                                                   experiment_desc
                                      organism_taxid  organism_name           library_strategy library_source  library_selection library_layout sample_accession sample_title instrument           total_spots total_size  run_accession run_total_spots run_total_bases
SRP178555       SRX5243234             GSM3560726: 97-gem-Day3-OEA115; Caenorhabditis elegans; RNA-Seq    GSM3560726: 97-gem-Day3-OEA115; Caenorhabditis elegans; RNA-Seq  6239            Caenorhabditis elegans  RNA-Seq          TRANSCRIPTOMIC  cDNA              SINGLE         SRS4245909                    Illumina HiSeq 2500  19052916     446793895  SRR8435699    19052916         971698716
SRP178555       SRX5243233             GSM3560725: 97-gem-Day3-OEA087; Caenorhabditis elegans; RNA-Seq    GSM3560725: 97-gem-Day3-OEA087; Caenorhabditis elegans; RNA-Seq  6239            Caenorhabditis elegans  RNA-Seq          TRANSCRIPTOMIC  cDNA              SINGLE         SRS4245908                    Illumina HiSeq 2500  18156224     428457264  SRR8435698    18156224         925967424
.. [truncated]

```
> 1 votes

---
Tags: phylogenetics, e-utilities
---


title: Installing HiC-Pro issues

text:
# Question
Title: Installing HiC-Pro issues
I am setting up a HiChIPseq pipeline, but installing software necessary for the HiC side of things has been unbecoming.

Both my PI and I have tried to install Hi-C Pro onto our Mac Pro that we use as our compute server (though in a couple of weeks we hope to migrate to an EPYC 64 Core computer server of our own!).

There are three approaches that I have taken to install Hi-C Pro.

1. From Source

To install from source, you need to supply the paths to its dependencies. Believing we did a good job, we run it and find problems accessing bx-python.

No matter how many times we update pip and try to install bx-python and specify that we want it installed for python2.

2. From singularity

This is the one I tried. From what I understand this creates a VM with the necessary dependencies installed, and allows you to use the program without building. However, it appears to abort when conda tries to install numpy.

```
+ conda install -y -c anaconda numpy=1.16.3
Collecting package metadata (current_repodata.json): done
Solving environment: failed with initial frozen solve. Retrying with flexible solve.
Collecting package metadata (repodata.json): / Killed
ABORT: Aborting with RETVAL=255
Cleaning up...

```

I have seen the failed solving environment bug many times, but have no idea how to solve it. I have tried updating conda, but that doesn't solve it.

3. Install Someone's Premade Conda Environment

Running conda.sh give me this:

```
Solving environment: failed

ResolvePackageNotFound:
<Lots of packages...>

```

Sorry if this is sloppy.

# Answer
> 1 votes
If you're reading this, I can only assume that's because you had trouble installing HiC-Pro too. I was finally able to install from source and run the software, so here is how I did it, giving detail to the error that I was hit with and the clues I used to solve the puzzles and get through installing this program:

> From Source To install from source, you need to supply the paths to its dependencies. Believing we did a good job, we run it and find problems accessing bx-python.
>
> No matter how many times we update pip and try to install bx-python and specify that we want it installed for p

# Setup

I downloaded from source once again using `git clone`.

The I `cat`ed the config file to find the programs I needed to find the path of, which I did using `which`. I ensured that I was using python2.

I also used `pip2` to install pysam, numpy, and scipy. I needed to use sudo to install pysam, though. I also installed something called bx-python.

# Make Configure

Here is the error I get when I run `make configure`

```
Checking Python libraries ...
Traceback (most recent call last):
File "../scripts/install/check_pythonlib.py", line 24, in <module>
      import bx.intervals
  ModuleNotFoundError: No module named 'bx'

  Traceback (most recent call last):
    File "../scripts/install/check_pythonlib.py", line 28, in <module>
      raise ImportError('Error - bx-python cannot be imported')
ImportError: Error - bx-python cannot be imported
Can not proceed without the required Python libraries, please install them and re-run
make[1]: *** [configure] Error 1
make: *** [configure] Error 2

```

Looking at this error for the fifth time, I realize that the first part, and not the last part, is the important one. It's yelling at me that I don't have the bx module installed. It also gives me the path to the file that threw the exception, so I read it and saw that this was what was being run:

```
## bx
try:
 import bx.intervals
 if vcmp(bx.__version__, '0.5.0') == -1:
   raise ValueError('bx-python '+bx.__version__+' detected. Version >= 0.5.0 required')
 except ImportError:
   raise ImportError('Error - bx-python cannot be imported')

```

I notice that there's nothing fancy, the thing just needs to be imported, and its actual name is bx. Wonderful, I pip2 it. But... I hit another error in pip:

```
pyrsistent requires Python '>=3.5' but the running Python is 2.7.16

```

Apparently, pyrsistent dropped support for python2. Wonderful. I get around it by:

```
sudo pip2 install awsebcli pyrsistent==0.16.0

```

In the middle of this installation, Im reading the output and see:

`bx-python 0.8.9 has requirement six>=1.13.0, but you'll have six 1.11.0 which is incompatible.`

So:

`sudo pip2 install six==1.13.0`

At this point, I, without any more problems, `sudo pip2 install bx`.

I run `make configure` again, and I run into...

```
Checking Python libraries ...
Traceback (most recent call last):
File "../scripts/install/check_pythonlib.py", line 24, in <module>
      import bx.intervals
  ModuleNotFoundError: No module named 'bx'

  Traceback (most recent call last):
    File "../scripts/install/check_pythonlib.py", line 28, in <module>
      raise ImportError('Error - bx-python cannot be imported')
ImportError: Error - bx-python cannot be imported
Can not proceed without the required Python libraries, please install them and re-run
make[1]: *** [configure] Error 1
make: *** [configure] Error 2

```

huh? Didn't I fix this?? Something seems fishy. For whatever reason, my mind then wanders to the Makefile to see if anything there pops out, cause... what else would be calling that python file? In the main Makefile I see:

```
make -f ./scripts/install/Makefile CONFIG_SYS=$(CONFIG_SYS)

```

So I head on over to ./scripts/install/Makefile and see:

```
./scripts/install/install_dependencies.sh -c $(CONFIG_SYS) -o  $(realpath $(PREFIX))/HiC-Proo    _$(VNUM) -q

```

Nothing like a goose chase. So I head on over to install\_dependencies.sh, where I see:

```
# python
which python > /dev/null;
if [ $? != "0" ]; then
     echo -e "$RED""Can not proceed without Python, please install and re-run""$NORMAL"
     exit 1;
else
     pver=`python --version 2>&1 | cut -d" " -f2` #<- THIS LINE
     vercomp $pver "2.7.0"
     if [[ $? == 2 ]]; then
         echo -e "$RED""Python v2.7.0 or higher is needed [$pver detected].""$NORMAL"
         exit 1;
     fi
fi

```

I check that lines output, AND IT'S PYTHON3!! For some reason, it didn't give me that error message, it just kept going from that line.

So, I changed all instances of `python` with `python2` in that file.

I reran `make configure` and got `The required Python libraries appear to be already installed.` Bliss.

# Make Install

Next step was `make install`, and guess what? A roadblock:

```
Traceback (most recent call last):
File "setup.py", line 2, in <module>
from numpy.distutils.core import setup
  ModuleNotFoundError: No module named 'numpy'
make: *** [iced] Error 1

```

Hmm... well, it's probably the same thing. python3 being called where python2 should be. So I double checked the main Makefile and changed the setup.py call to use python2.

I run `make install` and it runs for awhile, chugging along using the g++ compiler. Then I hit an odd error:

`cp: /Users/bosslab/Downloads/HiC-Pro-master/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4/HiC-Pro_2.11.4: name too long (not copied)`

I'm guessing that the Makefile has a recursive path somewhere. I check the appropriate section and find:

`cp -Ri $(MK_PATH) $(INSTALL_PATH)`

I simply replace that line with the actual paths, with the install\_path being adjacent instead of inside the mk\_path.

Then I get hit with:

`HiC-Pro installed in /Users/bosslab/Downloads/HiC-Pro-master/HiC-Pro_2.11.4 !`

So running `$install_dir/HiC-Pro --help` returns the help. YAY!

# Leftovers

I installed it, but how about when I actually run it?

`Error: The 'ice' command is not in your path. Please check where the 'iced' python package has been installed and update your PATH !`

Ok:

```
$ sudo pip2 install iced
Requirement already satisfied: iced in /Users/bosslab/Library/Python/2.7/lib/python/site-packages (0.4.2)

```

...so, I set up a test to see whether it actually is installed:

```
$ python2
Python 2.7.16 (default, Nov  9 2019, 05:55:08)
[GCC 4.2.1 Compatible Apple LLVM 11.0.0 (clang-1100.0.32.4) (-macos10.15-objc-s on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> iimport iced
>>> quit()

```

It's already installed. But then I had an epiphany:

```
$ python3
Python 3.8.5 (default, Jul 21 2020, 10:48:26)
[Clang 11.0.3 (clang-1103.0.32.62)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> ipomport iced
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
ModuleNotFoundError: No module named 'iced'
>>> quit()

```

Aha! Somewhere (i don't know where though) python3 is being called. So I just installed that module for python3 with `sudo pip3 install iced`.

Then I had some trouble with pointing to the right bowtie index (just give the directory where the index exists) and then I hit my next error dealing with HiC-Pro itself:

```
Pairing of R1 and R2 tags ...
Logs: logs/dixon_2M/mergeSAM.log
make: *** [bowtie_pairing] Error 1

```

Checking the log file named in the error message:

```
/usr/local/bin/python /Users/bosslab/Downloads/HiC-Pro-master/HiC-Pro_2.11.4/scripts/mergeSAM. py -q 0 -t -v -f bowtie_results/bwt2/dixon_2M/SRR400264_00_R1_hg19.bwt2merged.bam -r bowtie_re sults/bwt2/dixon_2M/SRR400264_00_R2_hg19.bwt2merged.bam -o bowtie_results/bwt2/dixon_2M/SRR400 264_00_hg19.bwt2pairs.bam
   File "/Users/bosslab/Downloads/HiC-Pro-master/HiC-Pro_2.11.4/scripts/mergeSAM.py", line 26
   print "Usage : python mergeSAM.py"
    ^
SyntaxError: Missing parentheses in call to 'print'. Did you mean print("Usage : python mergeS AM.py")?

```

Missing parentheses? Oh my god, it's using python3 for some reason again. I grep the installation directory for a file that calls merge\_sam.py, and find config-system.txt. I go in there and just need to change the python path from `/usr/bin/local/` where python3 is installed to `/usr/bin/` where python2 is installed.

Lastly, I get an error during the ICE normalization step that numpy is not installed. Using my big fat brain, I isntall numpy to my python3 and call it a day. Now it all works as intended.

# Answer
> 1 votes
I could install HiC-Pro using conda when creating a separate environment (which I always do for side-projects), and then pull the conda build of HiC-Pro:

```
conda create --name HiCPro python=2.7
conda activate HiCPro
conda install -c davebx hicpro

```

Not sure what `conda.sh` is that you describe, sometimes tricky to get random code from GitHub running. Give it a go with a standard conda install as I described.

---
Tags: chip-seq, software-installation, hi-c
---


# Question
Title: Find intersection among multiple columns and find the intersected elements
file input file

```
a <- read.csv("Studies_studies.tsv",sep = "\t")

#https://stackoverflow.com/questions/43658982/cross-comparison-of-columns-of-the-same-data-frame
bb <- tcrossprod(table(melt(as.matrix(a))[-1])) * !diag(37)

dd<- setNames(cbind.data.frame(t(combn(names(a), 2)),
                               combn(names(a), 2, function(x) length(intersect(a[, x[1]],a[, x[2]])))),
              c("col1", "col2", "count"))

```

So far my code.

I can get the count between each pairwise intersection happening between the columns but Im not able to get the elements that is common or intersected between two columns.

**For example**
This `C0026998 C1879321 1310` I do find that there are 1310 elements common between those to ID but I'm not able to find the name of the ID that are column.

So that way if I extrapolate my problem what i would like to see is the pairwise intersection among all the elements and i would like to store the output dataframe.

Any suggestion or help would be really appreciated.

input file dput

```
structure(list(C0026998 = c(32621177L, 32586363L, 32504186L,
32391628L, 32333156L, 32171069L), C1879321 = c(27276561L, 29286103L,
29988143L, 30076173L, 30024784L, 30892988L), C0023470 = c(27276561L,
29988143L, 30076173L, 30024784L, 26837842L, 29286103L), C0280449 = c(30024784L,
30076173L, 30892988L, 31688884L, 30651561L, 31516032L), C0522631 = c(26055299L,
12801839L, 10516750L, 28825596L, 12652469L, 16203817L), C1292775 = c(28883080L,
10400416L, 22581002L, 1503929L, 9499663L, NA), C0023472 = c(26837842L,
29091516L, 29498925L, 20525995L, 27217448L, 20525993L), C1292773 = c(30135184L,
24440648L, 26273061L, 21508125L, NA, NA), C4545381 = c(30093401L,
29933073L, 29437791L, 31227358L, 30541745L, 28555084L), C4707228 = c(26492932L,
21791467L, 28179273L, 24077845L, 30287390L, 22362038L), C2939461 = c(28588020L,
29988143L, 30076173L, 30510081L, 30024784L, 30963592L), C1336735 = c(28645776L,
28844816L, 31400961L, 30545576L, 29241450L, 28588020L), C2827362 = c(29119847L,
29486661L, 29567772L, 31775485L, 30670826L, NA), C3839741 = c(28494506L,
26492932L, 30923103L, 25715404L, 28556489L, 25006131L), C3839868 = c(28473620L,
25512507L, 31015209L, 29933069L, 29321554L, 21791467L), C4520840 = c(28246192L,
23935018L, 21750082L, 10741700L, 9407724L, NA), C4528668 = c(28645776L,
28588020L, 30510081L, 30409776L, 30093401L, 28844816L), C2825139 = c(29241450L,
30541745L, 31152020L, 24663049L, NA, NA), C1332153 = c(24440648L,
26273061L, NA, NA, NA, NA), C4721505 = c(18623376L, 21618422L,
21685468L, 25196856L, 16985182L, 18353632L), C2826177 = c(27276561L,
29286103L, 25512507L, 29079128L, 28473620L, 26789727L), C4324477 = c(29286103L,
22375971L, 30153096L, 30368038L, 25006131L, 30054307L), C2826176 = c(27276561L,
29988143L, 27959731L, 29286103L, 25512507L, 28588019L), C2827356 = c(27716285L,
29177434L, 30952642L, 30418178L, 27126994L, 29656442L), C3711382 = c(30510081L,
26492932L, 29991558L, 27927766L, 28490572L, 29860938L), C4704767 = c(30555165L,
30055822L, 25540937L, 30135184L, 30545923L, 30361063L), C0027013 = c(22375971L,
28962635L, 22375970L, 29515114L, 30185431L, 23216616L), C0153886 = c(30510081L,
28644114L, 28588020L, 30555165L, 30024784L, 29304833L), C0279094 = c(28644114L,
28588020L, 30510081L, 30555165L, 30076173L, 30024784L), C0279623 = c(27276561L,
30555165L, 28588020L, 30076173L, 30024784L, 30510081L), C1142169 = c(28588020L,
30510081L, 25512507L, 25540937L, 30409776L, 30668890L), C1275661 = c(29660836L,
30963592L, 22579233L, 21880637L, 26772158L, 27121471L), C1292776 = c(30055822L,
30115541L, 10025899L, 31152020L, 26577691L, 23129738L), C1320638 = c(27959731L,
18489989L, 25425683L, 28179273L, 30510079L, 30655375L), C1737261 = c(29988143L,
28645776L, 29933073L, 30545923L, 30361063L, 28336527L), C2826178 = c(25512507L,
29286103L, 28473620L, 25611491L, 21242187L, 19059939L), C4744720 = c(30555165L,
28644114L, 29597002L, 30541745L, 30118897L, 29660836L), relapse = c(32786187L,
32601798L, 32564196L, 32391628L, 32333156L, 32062741L)), row.names = c(NA,
6L), class = "data.frame")

```

output what i have

```
dput(head(dd))
structure(list(col1 = structure(c(3L, 3L, 3L, 3L, 3L, 3L), .Label = c("C0023470",
"C0023472", "C0026998", "C0027013", "C0153886", "C0279094", "C0279623",
"C0280449", "C0522631", "C1142169", "C1275661", "C1292773", "C1292775",
"C1292776", "C1320638", "C1332153", "C1336735", "C1737261", "C1879321",
"C2825139", "C2826176", "C2826177", "C2826178", "C2827356", "C2827362",
"C2939461", "C3711382", "C3839741", "C3839868", "C4324477", "C4520840",
"C4528668", "C4545381", "C4704767", "C4707228", "C4721505", "C4744720"
), class = "factor"), col2 = structure(c(18L, 1L, 7L, 8L, 12L,
2L), .Label = c("C0023470", "C0023472", "C0027013", "C0153886",
"C0279094", "C0279623", "C0280449", "C0522631", "C1142169", "C1275661",
"C1292773", "C1292775", "C1292776", "C1320638", "C1332153", "C1336735",
"C1737261", "C1879321", "C2825139", "C2826176", "C2826177", "C2826178",
"C2827356", "C2827362", "C2939461", "C3711382", "C3839741", "C3839868",
"C4324477", "C4520840", "C4528668", "C4545381", "C4704767", "C4707228",
"C4721505", "C4744720", "relapse"), class = "factor"), count = structure(c(1310L,
1310L, 164L, 9L, 6L, 23L), .Dim = 6L)), row.names = c(NA, 6L), class = "data.frame")

```

from the above output i do know that `C0026998 C1879321 1310` between C0026998 and C1879321 there are 1310 elements are common. Now the next thing i would like to find out is what are the elements that are common between those two IDs same goes for rest of the column which are to be compared pairwise.

The output I'm expecting is this

```
structure(list(`intersect(a$C0026998, a$C1879321)` = c(32621177L,
32586363L, 32504186L, 32391628L, 32333156L, 32171069L)), row.names = c(NA,
6L), class = "data.frame")

structure(list(`intersect(a$C0026998, a$C0280449)` = c(32171069L,
31841594L, 31770437L, 31753773L, 31688884L, 31516032L)), row.names = c(NA,
6L), class = "data.frame")

```

So in each column contain one Mesh term where there PMIDS which i want to find if these studies are overlapping or not.So in the above example i have given intersection between first four columns so i can see the IDs which are common. So i want to do the same with all the columns and save the results.

# Answer
> 2 votes
A tidyverse answer.

```
library("tidyverse")

results <- df %>%
  {gtools::combinations(ncol(.), 2, colnames(.))} %>%
  as_tibble(.name_repair="unique") %>%
  rename(col1=1, col2=2) %>%
  rowwise %>%
  mutate(intersection=list(intersect({{df}}[[col1]], {{df}}[[col2]]))) %>%
  unnest_wider(col=intersection, names_repair="unique") %>%
  pivot_longer(starts_with("..."), names_to="col", values_to="val") %>%
  pivot_wider(names_from=c(col1, col2), names_sep="_", values_from=val) %>%
  select(!col)

```

The first few results.

```
> results[1:5, 1:5]
# A tibble: 5 x 5
  C0023470_C00234~ C0023470_C00269~ C0023470_C00270~ C0023470_C01538~
             <int>            <int>            <int>            <int>
1         26837842               NA               NA         30024784
2               NA               NA               NA               NA
3               NA               NA               NA               NA
4               NA               NA               NA               NA
5               NA               NA               NA               NA
# ... with 1 more variable: C0023470_C0279094 <int>

```

---
Tags: r
---


# Question
Title: Choosing Fisher's Exact or Binomial test for overrepresentation in PANTHER
The PANTHER website offers a tool to obtain the GO-based overrepresentation of a gene list (the analyzed list) versus a reference gene list. After entering both these lists, the tool asks to choose between the Fisher's Exact or Binomial tests.

I guess that the convenience for one or another test would depend on the number of genes to analyze. If that is the case, which would be the ideal or advisable rank of gene numbers in both analyzed and reference lists (or their ratio) for each type of test?

(In particular I am working with a list of around 150 versus a reference list of around 2000 genes).

# Answer
> 2 votes
For a guide on using the tests in PANTHER, see https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6519457/.

In answer to this specific question, the defaults are the best choice (Fisher's Exact Test, with FDR correction). Fisher's Exact test gives more accurate statistics for relatively small counts (e.g. when only a few genes in your list are associated with a particular GO term). The FDR is a better correction to use for these types of analysis in general, especially given the large overlap between classes.

The other options are there mostly to support replication of older analyses, e.g. to replicate a published analysis that used the PANTHER over-representation tool prior to 2018, when the binomial with Bonferroni correction was the default.

---
Tags: gene, statistics, go
---