cannin/biostars_qa huggingface dataset:
Question: Hi guys,
I have a list of genes and I would like to know if among my genes there are targets of known or at least clinically validated drugs in the sense that at least the drugs have to be in some phase of clinical validation. No computational predictions, no hypothesis but only clinical validity. Drugs validated in vivo will be ok. Do you have an idea about the available tools to achieve my goal?

Thank you in advance

Answer: You could have a look at the Therapeutic target database: https://db.idrblab.org/ttd/. I've used it in the past without any filter, but I've seen it is possible to search for target AND disease.
Question: Dear all,

I have one problem, it is only one condition for my easy script for counting of GC content per row

awk 'NR>1{n=length($1); gc=gsub("[gcGC]", "", $1); print gc/n}' $i

How to count length of row without N character.

for example:

Input: ACAGCTTGCNNNN => length= 9 Gc content=5/9

format of output is not important, only how to count it.

Thanks a lot

Answer: Normalized GC content per row sans N

awk '{gsub("N","");t=length();gsub(/[GC]/,"");print int((t-length())/t*100)/100}'

or, assumes sequence symbols are strictly ACTGN.

awk '{gsub("N","");t=length();gsub(/[AT]/,"");print int(length()/t*100)/100}'
Question: This is a basic question. I want to rearrange a table like this:

```
A 1
B 2
B 3
B 4
C 5
C 6
```

to:

```
A 1
B 2,3,4
C 5,6
```

Could anyone gives a hint. Either with R or python?

Answer: You can try something like this:

require(plyr)

dat = read.csv(textConnection(
"A,1
B,2
B,3
B,4
C,5
C,6"),header=F)

print_V2=function(df){
tmp=c(as.character(df[1,1]),rep(NA,9))
tmp[2:(length(df[,2])+1)]=df[,2]
tmp
}

ddply(dat,.(V1),print_V2)

> ddply(dat,.(V1),print_V2)
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
1 A 1
2 B 2 3 4
3 C 5 6

Yes, you can `collapse` that too, as brentp suggested:

print_V2=function(df){
c(as.character(df[1,1]),paste(df$V2, collapse=", "))
}

ddply(dat,.(V1),print_V2)

> ddply(dat,.(V1),print_V2)
V1 V2
1 A 1
2 B 2, 3, 4
3 C 5, 6
Question: Hi all

I'm wondering if there is a way to exclude regions in samtools view.

Now I awk and sed lines away in post processing via piping.

Wouldn't it save time and space to just pileup more than one region at once into a file?

Maybe this is a dumb question but I don't have the software now while I'm wondering so I can't experiment.

For example:

samtools view chrY:2709520-3709520, chrY:17443438-27443438 ǀ samtools mpileup...
samtools view chrY:2709520-3709520 AND chrY:17443438-27443438 ǀ samtools mpileup.
samtools view chrY:2709520-3709520: chrY:17443438-27443438 ǀ samtools mpileup...

...

Answer: samtools view takes an argument to only output the alignments overlapping the regions of a BED file:

-L FILE output alignments overlapping the input BED FILE [null]
Question: Dear all,

I've been using GNU parallel for a while and it works quite well.
However, I recently needed to run some very long commands and parallel complained that the command line was too long:

parallel: Error: Command line too long (223235 >= 131049) at input 0: cat /tmp/10110507/adadev34gv_213...

It is a bit weird since my shell seems supports commands longer than 131049:

$ getconf ARG_MAX
2621440

and had no trouble running a very long command:

$ perl -e 'print("true "."x"x10000000);' | bash
$ echo $?
0

Does anyone know why it is so low and/or how to change it?
thanks,

Answer: I've been thinking a bit more about this and I think I've found a way around it using `extglob`.

Right now I have the commands like:

cat /long/path/file1 /long/path/file2 /long/path/file3 /long/path/file4 /long/path/file5 [...] /long/path/file10000 | mafft
cat /long/path/file10001 /long/path/file10002 /long/path/file10003 /long/path/file10004 /long/path/file10005 [...] /long/path/file20000 | mafft

But if I use `extglob` I could have them much shorter:

cat /long/path/file?(1|2|3|4|5 ... |10000) | mafft
cat /long/path/file?(10001|10002|10003|10004|10005 ... |20000) | mafft

So I tried running it through `GNU parallel`:

shopt -s extglob
. `which env_parallel.bash`
cat file_with_commands.sh | env_parallel [...]

and it seems to work. :)

thanks all for your help...
Question: Hi,

I‘ve got a list of genomic regions in python, say:

list=["chr1: 100-300", "chr1:150-350", "chr1:200-400", "chr1:500-700", "chr1:600-800", "chr1:900-1000"]

Some of them overlap each other. I wanna create a new list contain subsets of genomic regions and without overlapping regions in each subset. The result should be something like:

new_list=[["chr1:100-300", "chr1:500-700", "chr1:900-1000"], ["chr1:150-350", "chr1:600-800"], ["chr1:200-400"]]

I created my own function but it tooks 20 mins to achieve the result, if the origin list contain 1000 regions.

Is there any quick way to do that? Please help.

Answer: I would build a recursive function

def disjoint(l, res=[]):
sl = []
n = 1
while n < len(l):
p,q = l[n-1:n+1]
if p[2] >= q[1] and p[0]==q[0]:
l.remove(q)
sl.append(q)
else:
n+=1
res.append(l)
if sl: disjoint(sl, res)
return res

Running it on your example:

l = [('chr1',100,300), ('chr1', 150, 350), ('chr1', 200, 400), ('chr1', 500, 700), ('chr1',600, 800), ('chr1', 900, 1000)]
print disjoint(l)

Gives you what you want:

[[('chr1', 100, 300), ('chr1', 500, 700), ('chr1', 900, 1000)], [('chr1', 150, 350), ('chr1', 600, 800)], [('chr1', 200, 400)]]
Question: Hi,

I getting the following error. When I looked at my fastq files, I found that reads in R1 and R2 files are in the same order. What is the cause of this error and how to fix it? Thanks a lot.

```
abyss-map: invalid option -- l
abyss-map: invalid option -- 6
abyss-map: invalid option -- 4
Try `abyss-map --help' for more information.
Mateless 0
Unaligned 0
Singleton 0
FR 0
RF 0
FF 0
Different 0
Total 0
abyss-fixmate: error: All reads are mateless. This can happen when first and second read IDs do not match.
error: `pe150_1-3.hist': No such file or directory
make: *** [pe150_1-3.dist] Error 1
make: *** Deleting file `pe150_1-3.dist'

real 0m1.379s
user 0m0.015s
sys 0m0.082s
```

Answer: This error seems to be a glitch in abyss v.1.3.6. When I used the latest abyss version (1.9.0), abyss-pe successfully submitted abyss-map jobs. So the solution to this problem is to use the latest version of abyss.
Question: The peaks of a particular Transcription factor CHIP-seq as identified by multiple peak calling algorithms e.g. MACS2 peakzilla e.t.c identify peaks of varying sizes or regions from 50 bp up to 1000 bp or more. In that case I was wondering that if the average size of transcription factor binding site is 10 bp based on available data in OReganno e.t.c (as indicated by [this post][1]) would it be safe to assume that for a particular peak the TFBS location is in its center or not. If not then what would be a good way to go about doing it and if such a method already exist or if someone has already done it.

[1]: https://www.biostars.org/p/64854/

Answer: You could use FIMO in the MEME suite to scan for motif models (JASPAR, etc.) across your genome of interest, at a level of statistical significance that you deem acceptable (*e.g.*, 1e-4, 1e-5).

Take the search result from FIMO and convert it to a UCSC BED-formatted file. This file contains the putative binding sites of all motifs from the motif models across the whole genome.

Then do set operations with BEDOPS tools (like *bedmap*) to precisely locate putative TF binding sites that overlap with - are contained entirely or partially within - your ChIP-seq peaks.

$ bedmap --echo --echo-map peaks.bed wg-motifs.bed > answer.bed

If you repeat your experiments with other motifs, you can reuse your whole-genome search result to apply the same set operations.
Question: What is the difference between paired end reads and overlapping reads, after reading around I am struggling to find an answer as to what overlapping reads actually are. From what I can see they just seem to be another term for paired end reads but I feel this must be wrong or I would have found somewhere stating this?

Secondly, why merge overlapping reads prior to assembly?

Any help would be appriciated :)

Answer: well, overlapping reads are a "special case" of paired end reads.

Just as the word say they overlap each other (while in normal cases that should not be). It's actually all linked to the fragment/insert size of your library and the read length sequencing mode you execute.

eg. say your insert is 200bp and you sequence 150bp PE reads, those reads will overlap by 50-100 bases (== overlapping reads), if for the same library you do 75bp PE sequencing they will not overlap.

looking at the image above, if the inner distance becomes negative, the reads will overlap

In some cases (eg. Assembly, amplicon sequencing, ... ) it's advantageous to have overlapping reads as you will end up with longer sequences which you know form a single biological entity
Question: Hi, For Gene enrichment analysis, I used ClueGO as a Cytoscape plugin. At first, I updated ontologies.then, to enrichment analysis for my 40 Ensemble ids, I executed ClueGO in Cytoscape. The process will be continued but in the final stage for getting the GO network, my Cytoscape is closed suddenly. I have done this process but the results were similar.

My system and Cytscape Info:

Cytoscape Version: 3.7.2
ClueGO v2.5.6 + CluePedia v1.5.6
Java: 1.8.0_231 by Oracle Corporation Java
Home: c:\program files\java\jre1.8.0_231
OS: Windows 10 10.0 - amd64

I appreciated it if anybody guides me in this problem.

Best Regards,

Answer: If Cytoscape is closing, it means that something is crashing. There are two things you can try:

1. You could upgrade your Java to Java 11 (a much more current release) and your Cytoscape to 3.8, which is our latest. Note that Cytoscape works only with Java 11 (not Java 1.8, 12 or 13).
2. Alternatively, you could try to run Cytoscape from the command prompt and see what error is getting produced before it closes. That will help us to debug it.

-- scooter
Question: I have a dataset of >100 different samples. Samples are from different genotypes (e.g. X, Y, Z) and 4 different time points (T0,1,2,3) with 3 biological replicates (R1,2,3). I'm measuring values for 50 different genes (in raws)

structure(list(Gene = structure(1:2, .Label = c("A", "B"), class = "factor"), X_T0_R1 = c(1.46559502, 0.220140568), X_T0_R2 = c(1.087642983, 0.237500819), X_T0_R3 = c(1.424945196, 0.21066267), X_T1_R1 = c(1.289943948, 0.207778662), X_T1_R2 = c(1.376535013, 0.488774258), X_T1_R3 = c(1.833390311, 0.182798731), X_T2_R1 = c(1.450753714, 0.247576125), X_T2_R2 = c(1.3094609, 0.390028842), X_T2_R3 = c(0.5953716, 1.007079177), X_T3_R1 = c(0.7906009, 0.730242116), X_T3_R2 = c(1.215333041, 1.012914813), X_T3_R3 = c(1.069312467, 0.780421013), Y_T0_R1 = c(0.053317766, 3.316414959), Y_T0_R2 = c(0.506623748, 3.599442788), Y_T0_R3 = c(0.713670106, 2.516735845), Y_T1_R1 = c(0.740998252, 1.444496448), Y_T1_R2 = c(0.648231834, 0.097957459), Y_T1_R3 = c(0.780499252, 0.187840968), Y_T2_R1 = c(0.35344654, 1.190274584), Y_T2_R2 = c(0.220223951, 1.367784148), Y_T2_R3 = c(0.432856978, 1.403057729), Y_T3_R1 = c(0.234963735, 1.232129062), Y_T3_R2 = c(0.353770497, 0.885122768), Y_T3_R3 = c(0.396091395, 1.333921747), Z_T0_R1 = c(0.398000559, 1.286528398), Z_T0_R2 = c(0.384759325, 1.122251177), Z_T0_R3 = c(1.582230097, 0.697419716), Z_T1_R1 = c(1.136843842, 0.804552001), Z_T1_R2 = c(1.275683837, 1.227821594), Z_T1_R3 = c(0.963349308, 0.968589683), Z_T2_R1 = c(3.765036263, 0.477443352), Z_T2_R2 = c(1.901023385, 0.832736132), Z_T2_R3 = c(1.407713024, 0.911920317), Z_T3_R1 = c(0.988333629, 1.095130142), Z_T3_R2 = c(0.618606729, 0.497458337), Z_T3_R3 = c(0.429823986, 0.471389536)), .Names = c("Gene", "X_T0_R1", "X_T0_R2", "X_T0_R3", "X_T1_R1", "X_T1_R2", "X_T1_R3", "X_T2_R1", "X_T2_R2", "X_T2_R3", "X_T3_R1", "X_T3_R2", "X_T3_R3", "Y_T0_R1", "Y_T0_R2", "Y_T0_R3", "Y_T1_R1", "Y_T1_R2", "Y_T1_R3", "Y_T2_R1", "Y_T2_R2", "Y_T2_R3", "Y_T3_R1", "Y_T3_R2", "Y_T3_R3", "Z_T0_R1", "Z_T0_R2", "Z_T0_R3", "Z_T1_R1", "Z_T1_R2", "Z_T1_R3", "Z_T2_R1", "Z_T2_R2", "Z_T2_R3", "Z_T3_R1", "Z_T3_R2", "Z_T3_R3"), class = "data.frame", row.names = c(NA, -2L))

I want to average reach 3 replicates (R1,2,3) for each time point (T0,1,3,5) and make a new matrix and then create a heat map with dendrogram.

How can I average each 3 replicates of each genotype at particular time point and make the new matrix? is this possible with tydr?

Answer: since no one has replied, following is the answer gene vs time heatmap. Data frame (df) is constructed from OP data and replicates are average (mean).

library(tidyr)
library(dplyr)
library(ggplot2)

df2=gather(df,k,expression,-Gene) %>%
separate(k,c("sample","time","replicate"))%>%
group_by(Gene,sample,time) %>%
summarise(mean=mean(expression))

ggplot(data=df2,aes(time, Gene, fill=mean)) +
geom_tile()+
scale_fill_gradient(low = "green", high = "red")
Question: Hi,

How can I get the Internal Stability graph or short oligonutides, like this one for B1 ([1993 Rychlik][1])

![enter image description here][2]

Also I don't know why the last nucleotides does not match with any value of the graph

I try to get the graph using nearest neighbor dG values, but I cant get a similar graph with the python 2.7 code below, using this approach I can get every value exept for the last nucleotide.

![enter image description here][3]

def calc_dG(seq):
NearPairsPlus = ['AA','AC','AG','AT',
'CA','CC','CG','CT',
'GA','GC','GG','GT',
'TA','TC','TG','TT']
DeltaG = [-1.9, -1.3, -1.6, -1.5,
-1.9, -3.1, -3.6, -1.6,
-1.6, -3.1, -3.1, -1.3,
-1.0, -1.6, -1.9, -1.9 ]
c = []
seq = seq+ ' '
a = 0
while a < len(seq):
b = 0
nuc = []
while b < len(NearPairsPlus):
if seq[a-2:a] == NearPairsPlus[b]:
c.append(DeltaG[b])
b = b + 1
a = a + 1
return c

sequence = "ACTTGGGATTGGGCT"
dG_values = calc_dG(sequence)
y = dG_values
x = list(sequence)[:-1]

import matplotlib.pyplot as plt

plt.plot(dG_values)
plt.gca().invert_yaxis()
plt.ylabel(r'$\Delta$ G $\left(\frac{kcal}{mol}\right)$')
plt.show()

[1]: https://www.ncbi.nlm.nih.gov/pubmed/21400260
[2]: https://i.snag.gy/DBjuKk.jpg
[3]: https://i.cubeupload.com/FC5rQv.png

Answer: I just solved it, didn't read that the dG values were from pentamers.

Now I got a similar graph from B1
![enter image description here][1]

def calc_dG(seq):
NearPairsPlus = ['AA','AC','AG','AT',
'CA','CC','CG','CT',
'GA','GC','GG','GT',
'TA','TC','TG','TT']
DeltaG = [-1.9, -1.3, -1.6, -1.5,
-1.9, -3.1, -3.6, -1.6,
-1.6, -3.1, -3.1, -1.3,
-1.0, -1.6, -1.9, -1.9 ]
c = []
seq = seq+ ' '
a = 0
while a < len(seq):
b = 0
nuc = []
while b < len(NearPairsPlus):
if seq[a-2:a] == NearPairsPlus[b]:
c.append(DeltaG[b])
b = b + 1
a = a + 1
return c

sequence = "ACTTGGGATTGGGCT"
dG_values = calc_dG(sequence)
y = dG_values
x = list(sequence)[:-1]

dG_pentamers = []
for i in range(0, len(y)):
dG_pentamers.append(sum(y[i:i+4]))
if i+4 == len(y):
break

import matplotlib.pyplot as plt

y_labels = list(sequence[0:-5])
y_labels.append(sequence[len(sequence)-5:len(sequence)])

plt.plot(dG_pentamers, marker='o', linestyle=':', color='b')

plt.ylabel(r'$\sum \Delta$ G $\left(\frac{kcal}{mol}\right)$ pentamers')
plt.xticks(range(len(y_labels)), y_labels)
plt.gca().axes.set_ylim([-12,-5])
plt.gca().invert_yaxis()
plt.show()

[1]: http://i.imgur.com/AhtXrT5.png
Question: Hi Everyone!!

**Background:** I got a unique single cell dataset produced in house, where there is low number of genes expressed (initial stage of developmental cycle). We have control and drug treatment scRNASeq samples and post differential expression analysis, I wish to perform GO-Term analysis. However, I am confused about how to create an appropriate background (Universe Gene-set) and an appropriate methodology to carry out GO Term analysis

**Method 1 for Background creation**

Find genes which are expressed even in 1 cell in both control and treatment and then take union of Gene IDs.

**Method 2 for Background creation**

Find genes which are expressed in either 10% of cells in Control or 10% of cells in treatment. For this I am trying to calculate proportion of cells expressing each gene in each sample i.e. in ctrl and treatment separately using this function:

```r

per.gene.per.sample.pct <- function(sobj,sample_col){
## Making the table summarising number of cells in each sample
i=sobj@meta.data %>% plyr::count(sample_col)
rname=rownames(sobj[["RNA"]]@counts)
## Calculating the percentage of cells where the gene is expressed
## in each sample

rowsum.exp <- function(x){
rowSums2(sobj[,sobj@meta.data[,sample_col]==i[,1][x]][["RNA"]]@counts>0)/i[,2][x]}
tt <- 1:nrow(i) %>% purrr::map(function(x) rowsum.exp(x))
tt<- t(plyr::ldply(tt))
rownames(tt) <- rname
colnames(tt) <- i[,1]
return(tt)
}

```

Besides, for GO Term analysis should I combine all Differential genes from all the clusters and then carry out GO term analysis or should I separately carry out GO Term analysis for each cluster.The aim is to find out eventually which pathways got perturbed due to treatment and not the cell types in single cell data.

**PS:** When I try Method 2 for background creation, It leaves me with very few genes as compared to number of DE genes that I found using Seurat!! I understand this since the DE genes are calculated on per cluster basis and every gene need not be necessarily expressed in every cluster. But then does it mean that I should make separate background geneset based on genes expressed in each cluster?

On the other hand, Method 1 gives me exorbitant number of genes and I think that's not ideal either.

Answer: I think method 1 should work just fine, but maybe you can set a cell number threshold. Eg must be expressed in at least 10 cells in a cell state. This way the number of detected genes should increase as more cells are added to the dataset. Let us know how many genes are detected with this approach and how many genes are in the background and foreground sets.

> Besides, for GO Term analysis should I combine all Differential genes
> from all the clusters and then carry out GO term analysis or should I
> separately carry out GO Term analysis for each cluster.The aim is to
> find out eventually which pathways got perturbed due to treatment and
> not the cell types in single cell data.

I think treating the different cell states as separate "experiments" is typical. We have to be open to the idea that the different cell states will respond differently to the disease/stimulus.
Question: I have my data as "10x" format.

```
barcodes.tsv
genes.tsv
matrix.mtx
```

How can I convert this to .cloupe format to view in the 10x loupe browser?

Answer: I asked 10x once...you can't. You need to run cellranger on the fastqs to make a loupe file.
Question: I have downloaded the file "protein chemical.links.v4.0.tsv.gz" from STITCH 4.0. From the database, the chemical identifiers are described as below:

**Chemical identifiers are of the form CID123456789, where the number stands for the PubChem compound identifier of the merged compound**

Here are some example of chemical identifiers in the file :

```
CID171549452
CID171549315
CID171464722
```

However, I can't find the result in PubChem database by inputting the identifiers I obtain from the downloaded file. Can I know what are the problems behind and how can I solve this?

Answer: The description from STITCH is not entirely accurate. STITCH has two kinds of chemical identifiers CID0s and CID1s. The CID0s are individual compounds, CID1s are the merged compounds. The number being referred to is the number after CID0 or CID1. CID171549452 is thus the merged compound represented by compound 71549452 in PubChem.

In the coming version of STITCH we will change the CID0/CID1 terminology with CIDs/CIDm, where `s` is single compounds and `m` is merged compounds.
Question: I just found out that TreeView is great to have a general heatmap of our RNAseq samples (Transcriptome). I installed TreeView on my computer and had a look at there sample files, though their website isn't really helping and I don't know where to start.

http://jtreeview.sourceforge.net/

**Is there any good manual or tutorial to get started with TreeView**?

Answer: Its straight forward with usage and also it has nice [manual][1]. Once you get your expression data (generally log transformed and normalized) you cluster them using [cluster 3.0][2]. Cluster generates mainly 3 files.

1. .cdt (cluster data) file which is the one you open with TreeView. It contains the same expression table that you used for clustering except the data has be re organised according to clustering.
2. .gtr and .atr (gene tree and array tree) contain structure of your expression data and your arrays (in this case RNA-seq replicates and samples).
3. If you are using k-mean clustering you will have .kgg file.

Once you are done with clustering, you open .cdt file with TreeView. If the heatmap which shows up is not visible clearly, go to settings -> pixel settings and change contrast.

You can search for your gene in analysis tab with find genes. You can export the image as png or jpeg. You can also do PCA and other stuffs, you will have to play around.

Remember that cluster generally crashes with low memory computer when you have thousands of gene.

[1]: http://jtreeview.sourceforge.net/docs/JTVUserManual/single.html
[2]: http://bonsai.hgc.jp/~mdehoon/software/cluster/software.htm
Question: I want to calculate **markerwise** Weir And Cockerham Fst for certain differentiating markers. I have my data in PLINK(bed/bim/fam) format. I have scourged the internet for packages to calculate the fst but very few seem to support native plink format directly as an input.

There exists **snpstats** (bioconductor R package) but the method of calculation seems to be a non standard one.

Can anyone suggest an R package/software that can calculate Fst(preferebly Weir And Cockerham) from **PLINK format data** as input?

Answer: How are your subpopulations defined? (I may just go ahead and add this to PLINK, but I first want to make sure that your subpopulations are defined in a PLINK-friendly manner.)
Question: Hi,

I used CCAT to call peaks for the H3k4me3 mark in my samples and the average peak width I noticed is ~3kb

I wanted to get a sense check on these nos. ? Usually in cancer samples, is that the sort of profile that is seen for H3k4me3. I tried looking this up in literature but I dint get enough data pts and the few that I did showed an average profile of 1.5kb.

Does anybody have any references on this ?

Thanks !

Answer: From H3K4me3 in mESC (mouse embryonic stem cells), what I know is that the promoters which have the tri-methylated nuclesomes at active. This width, I plot as an average profile in and around +/-3KB of TSS, which is sufficient to answer. Its hard to say the exact width, but I would assume it to be a max of 3KB so it might be fine.

As an example, I took this screenshot from the cistrome webpage. Its obvious that the average width would be around 2-3KB but I am not sure if that be different in the cancer cell line.