cannin/biostars_qa huggingface dataset:
Question: Hi,

I am trying to run a data of 105 genes (combined SNP,INDEL,CNV) in the cbioportal webpage to understand how these genes are affected across Ovarian serous cystadenocarcinoma data(Nature, 2011). First of all I cannot do it for more than 100 genes. So I split the list in 2. But most improtantly I see that I cannot download the oncoprint image in pdf. It seems it does not work. If I try to mess much with the page clicking different option it crashes most of the time. I can only download in SVG format as of now. I would like to ask that can the sorting be done based on the frequency percentage? If that is possible how can that be done? I f I try to drag the genes with high frequency percentage up to the top the page crashes. I have tried in Mozilla and Chrome. Pages are either crashing or becoming unresponsive. I do not see that option in the webpage. Any help appreciated

Regards

Answer: There is some technical issues for downloading large PDF files, please use SVG and convert to PDF using tools such as Adobe Illustrator.

There is no functions for sorting genes on cBioPortal. I would suggest to get the percentage first and then copy the results of alteration percentages per gene in the Download tab and sort them in a spreadsheet. Then you could do another query with genes sorted by alteration frequency.
Question: Hi everyone.

I found Different ENSG IDs for the same gene_name in gtf file (gencode.v26.annotation.gtf).

For example:
"ENSG00000273544.1" and "ENSG00000252840.1" both have a gene_name "SNORA44".

I searched for these ENSG IDs on ensemble, and got 2 answer:
ENSG00000273544.1 is a 132bp snoRNA at Chromosome 1: 28,580,381-28,580,512;
ENSG00000252840.1 is a 108bp snoRNA at Chromosome 1: 151,527,831-151,527,938;

So why do they share the same gene_name?

Thank you in Advance.

Answer: These cases happen very often. In this particular instance, notice that you are using an old version of GENCODE, 26.

If you check the latest [GENCODE 39][1], you can see that ENSG00000273544.1 is the correct entry for this gene symbol. You can also check [Gene Cards][2] where you will see that there is only one ensembl id. Finally, the [ensembl page][3] for the other id shows that ENSG00000252840 is no longer associated with the symbol SNORA44.

[1]: https://www.gencodegenes.org/human/
[2]: https://www.genecards.org/cgi-bin/carddisp.pl?gene=SNORA44
[3]: https://www.ensembl.org/Homo_sapiens/Gene/Summary?db=core;g=ENSG00000252840;r=1:151527831-151527938;t=ENST00000517031
Question: I am trying to do the example in [this link][1].

[1]: http://wgs-assembler.sourceforge.net/wiki/index.php/PBcR

the lambda example, I use this command:

```
wgs-8.2alpha/Linux-amd64/bin/PBcR -length 500 \
-partitions 200 \
-l lambdaIll \
-s pacbio.spec \
-fastq pacbio.filtered_subreads.fastq \
genomeSize=50000 \
illumina.frg
```

but it gives me this error:

```
Running with 40.01948X (for genome size 50000) of lambdaIll sequences (2000974 bp).
Correcting with 50X sequences (2500000 bp).
AMOS or PBDAGCON is required to call consensus and neither was available. Please add either SMRTportal or AMOS to your path and try again
```

I already have *pbdagcon* in my path.

Any Idea

Thanks in advance

**Update**: after installing AMOS It worked

Answer: Right now I am installing AMOS and its dependency and I will see what will happens and your suggestion is to use `-pbCNS` flag ?!

**Update**: After installing AMOS It worked
Question: After reading the paper and docs, I am having a little trouble understanding the difference and usage for calling at the segment level or the bin level. I am running CNVkit using human exome seq data.

My questions are:

What use case would the segment level calls be best for

What use case would the bin level calls be best for

Which one is more accurate in what context

What is the point of having two? Is the bin level output noisier but at a higher resolution than the segment level? Just assuming default segmentation and bin sizes. I'm really confused here -- anything helps!

Answer: The two file types might be more intuitive if you have experience with an older microarray-based method, [array comparative genomic hybridization][1] (aCGH). Bins are equivalent to microarray probes there.

The bins provide a fine-grained genome-wide copy number signal plus some noise. Segmentation attempts to remove the noise and infer the location of discrete copy number alterations, i.e. the individual regions that have been duplicated or deleted. Segments that are not neutral (i.e. diploid, log2=0) are putative copy number alterations. CNVkit's `call` command helps infer more about the segments beyond their breakpoints.

So, in general, use the segments (.cns) for most follow-up analysis. The bin-level data (.cnr) is useful for plotting and showing the level of support for each segment, and for tracking down potential artifacts like especially noisy regions of the genome. Also, in the case of small or single-exon CNA, the .cns file will typically not include it but the .cnr file may show some evidence for the copy number change that you could then look to confirm independently.

[1]: https://en.wikipedia.org/wiki/Comparative_genomic_hybridization
Question: Hello,
I have used DAVID tool for the Functional annotation of my data (microarray data) and now i wan to visualize the results of GO term in the form of Bar plots (i have seen the such results from the research paper)[reference image from the paper][1]
the link contains the image of box plot in the research paper.
My question is,
1. Which tools can be used to create such plots or any online tool
I have used Revigo tool but it doesn,t give such kind of plots.
Any suggestions would be helpful.

[1]: https://www.spandidos-publications.com/ijo/44/3/858

Thank you.

Answer: Hi, it is just a simple barchart or bar-plot. First, create a data-frame that has your GO enrichment terms as rows and count values in a single column, like this:

input
GO1 128
GO2 1
GO3 70
GO4 6
GO5 17
GO6 14
GO7 18
GO8 4
GO9 13
GO10 0
GO11 0
GO12 5
GO13 16
GO14 5
GO15 2
GO16 0
GO17 0
GO18 20
GO19 6
GO20 0
GO21 0
GO22 0
GO23 2
GO24 0
GO25 6
GO26 1
GO27 5
GO28 1
GO29 0
GO30 0
GO31 0
GO32 0
GO33 1

Then, generate separate barcharts with *lattice* and piece them together on the same panel using *grid* / *gridExtra*:

require(lattice)
require(RColorBrewer)
require(rasterVis)
require(gridExtra)
require(grid)

A <- barchart(
data.matrix(input[1:15,]),
col="black",
xlab="Total count", ylab="",
pretty=TRUE,
scales=list(
x=list(cex=1.0, col="black", fontface="bold"),
y=list(cex=0.8, labels=rownames(input), col="black", fontface="bold"),
xlab=list(cex=2.0),
ylab=list(cex=2.0),
tck=c(1,0)),
aspect="fill")

B <- barchart(
data.matrix(input2[1:15,]),
col="black",
xlab="Total count", ylab="",
pretty=TRUE,
scales=list(
x=list(cex=1.0, col="black", fontface="bold"),
y=list(cex=0.8, labels=rownames(input2), col="black", fontface="bold"),
xlab=list(cex=2.0),
ylab=list(cex=2.0),
tck=c(1,0)),
aspect="fill")

layout <- cbind(c(1,2))
grid.arrange(
layout_matrix=layout,
arrangeGrob(A,
top=textGrob("A",
x=unit(0.05,"npc"),
y=unit(0.9,"npc"),
just=c("left","top"),
gp=gpar(fontsize=32))),
arrangeGrob(B,
top=textGrob("B",
x=unit(0.05,"npc"),
y=unit(0.9,"npc"),
just=c("left","top"),
gp=gpar(fontsize=32))),
ncol=1)

------------------------------------------------------------
---------------------

For side to side, alter the layout and use ncol=2 to grid.arrange

A <- barchart(
data.matrix(input),
col="black",
xlab="Total count", ylab="",
pretty=TRUE,
scales=list(
x=list(cex=1.0, col="black", fontface="bold"),
y=list(cex=0.8, labels=rownames(input), col="black", fontface="bold"),
xlab=list(cex=2.0),
ylab=list(cex=2.0),
tck=c(1,0)),
aspect="fill")

B <- barchart(
data.matrix(input2),
col="black",
xlab="Total count",
ylab="",
pretty=TRUE,
scales=list(
x=list(cex=1.0, col="black", fontface="bold"),
y=list(cex=0.8, labels=rownames(input2), col="black", fontface="bold"),
xlab=list(cex=2.0),
ylab=list(cex=2.0),
tck=c(1,0)),
aspect="fill")

layout <- cbind(c(1),c(2))
grid.arrange(
layout_matrix=layout,
arrangeGrob(A,
top=textGrob("A",
x=unit(0.05,"npc"),
y=unit(0.9,"npc"),
just=c("left","top"),
gp=gpar(fontsize=32))),
arrangeGrob(B,
top=textGrob("B",
x=unit(0.05,"npc"),
y=unit(0.9,"npc"),
just=c("left","top"),
gp=gpar(fontsize=32))),
ncol=2)
Question: Hey everyone. I am currently working on automating a SNP pipeline and need some help figuring out how to interface Python with GenomeStudio. I want to be able to pass input from a python script into GenomeStudio so that the entire analysis process can be automated and manual intervention is only required rarely. Does anyone have any experience with this? I've done a few Google searches and haven't come up with anything.

Thanks,
Travis Couture

Answer: Not exactly automating GenomeStudio, but what you could do is turn on AutoConvert on the iScan instrument and have it generate gtc files for each sample in your experiment. Then you can use the recently published open-source Python gtc file parser from Illumina to write a script for generating reports or whatever else you want to do with the data: https://github.com/Illumina/BeadArrayFiles. In the examples folder of that repo, there's an example script that lets you generate final report files from gtc files, thereby letting you avoid having to use GenomeStudio at all.

Or if you're working with idat files instead, what you could do is install and run AutoConvert off the instrument and write some Python code to go from idat + bpm + egt --AutoConvert--> gtc --python gtc file parser--> reports
Question: I'm reading up on the Seurat user guide: https://satijalab.org/seurat/v3.1/pbmc3k_tutorial.html
And they mention for QC utilizing

**The number of unique genes detected in each cell** and **The total number of molecules detected within a cell**

They then refer to them as nCount_RNA and nFeature_RNA, but I'm not sure which is which. So my question is:

1.) What are the nCount_RNA and what are the nFeature_FNA
2.) Later in the pipeline, when you're normalizing the data, it says they "normalizes the feature expression measurements for each cell by the total expression." Can anybody explain that?

Answer: `nFeature_RNA` is **the number of genes detected in each cell**. `nCount_RNA` is **the total number of molecules detected within a cell**. Low `nFeature_RNA` for a cell indicates that it may be dead/dying or an empty droplet. High `nCount_RNA` and/or `nFeature_RNA` indicates that the "cell" may in fact be a doublet (or multiplet). In combination with %mitochondrial reads, removing outliers from these groups removes most doublets/dead cells/empty droplets, hence why filtering is a common pre-processing step.

The `NormalizeData` step is basically just ensuring expression values across cells are on a comparable scale. By default, it will divide counts for each gene by the total counts in the cell, multiply that value for each gene by the `scale.factor` (10,000 by default), and then natural log-transform them.
Question: Hello,

I am using HOMER to find enriched known motifs in my differentially expressed gene lists (RNA-seq). I am using the findMotifs.pl tool. I get enriched motifs, and I would like to know now which genes have this motif. Is it possible to get these genes? I can't find it in the manual pages of HOMER, did I miss it? Does anyone has an idea how to get the target genes?

Ben

Answer: Sorry for posting this question, I found what I wanted in the manual pages section:

Finding Instances of Specific Motifs
Question: Hello all. I would appreciate your help with this. Im doing an WCGNA, everything is getting well but the only thing that i have a trouble is when i tried to generate a file with all my genes from each module (to perform a gene ontology analysis)
Ive performed this command, it doesnt tell me there is an error, but when i open my .txt files those are all empty.
What would you could suggest for me?
Thank you so much

allLLIDs = pulmonData$gene_symbol
intModules = unique(moduleColorsAutomatic)
for (module in intModules)
{
modGenes = (moduleColorsAutomatic==module)
modLLIDs = allLLIDs[modGenes];
fileName = paste("LocusLinkIDs-", module, ".txt", sep="");
cat(modLLIDs, sep="\n", file=fileName)
}

TotfileName = paste("LocusLinkIDs-all.txt", sep="");
cat(allLLIDs, sep="\n", file=TotfileName)

Answer: ive already answered it.

The problem was that the `pulmonData` first column was named "Gene ID", i changed to "gene_symbol" and now it works.

Thank you!
Question: Hi everyone. I am dealing with a parasite that contains `~5000` genes. The scRNAseq data is from a stage of its life cycle where few genes are expressed (the median genes expressed are ~100 in our data).

1. In that case, does it make sense to use `nfeature` default value `2000` while running `SelectIntegrationFeatures `.
2. Also, the UMAP looks as given below (we have data from two-time points with 1 hr difference between T1 and T2). So should I perform integration or not?
3. Finally since there is no profound batch effect should I still use `SCTransform` or follow the `pbmc3K` workflow?

I am not getting many DE genes either way (with or without integration) and DE genes number fairly overlap.

![enter image description here][1]

[1]: /media/images/c3e82de0-d1ad-41e4-bcff-1cf6c6eb

Answer: I'd probably just try it with a few different numbers and see how it does (e.g. 2000, 1000, 500, 250, 100, etc). Hard to answer your second question without just slapping both timepoints in the same plot. I'd probably still try integration.
Question: I want to download a SRA data (SRA052842) with the prefetch command in command line of Ubuntu, The prefetch command is belonged to SAM tools package, so I've used the following command:

prefetch SRA052842

and I've faced with following error:

Maximum file size download limit is 20,971,520KB

2015-11-07T06:28:05 prefetch.2.3.5 err: path not found while resolving tree within virtual file system module - 'SRA052842' cannot be found.

Where is the problem? How can I fixed it?

Answer: You have two choices

1. Use http://www.sradownload.com

2. Use a FTP client program such as Filezilla under an anonymous identity
- Just download and install the program
- Create a new connection pointing at ftp-trace.ncbi.nhi.com as the server
- Hit the advanced tab
- Fill the remote directory with /sra/sra-instant/reads

And you will be allowed to choose your SRA data by Exp, Run, Sample and Study
Question: Hello Forum,

What pipelines are available for exome sequencing analysis? Thanks for information.

Answer: Probably in need of an update, but not a bad overview of the steps: http://seqanswers.com/wiki/How-to/exome_analysis
Question: Hi,
I have a set of piRNAs in fasta format and I want to find their potential targets.
According to literature nucleotides in position 1-10 in the piRNAs sequences should constitute the seed and do not tolerate any mismatch while nt 11-20 tolerate up to 3 mismatches.
Bowtie can specify the length of the seed, but I'm having troubles finding the right combination of parameters to translate the piRNA target pairing logic into a single bowtie command line.

Is there any solution to get this done with a single alignment run or do I have to run a first alignment guided by the seed 1-10nt, and successively another alignment for the 11-20 region allowing for up to 3 mismatches?

Thank you in advance for your help!

Answer: Or maybe you can use [LAST](http://last.cbrc.jp) with a custom [seeding scheme](http://last.cbrc.jp/doc/last-seeds.html) ?
Question: Trying to filter loom file in R.

I have done analysis with seurat and trying to add another layer with velocity analysis scvelo (python). I could probably also do it in python with loompy library.
From what I saw, the velocity analysis start from the 10x directory but I have already filtered seurat object and save it in an rds file. I have a list of ids in a csv files. I want only use these cells for the analysis of velocity.

I'm using ReadVelocity to read the loom file and then convert to seurat but then I'm trying to subset the object .

I'm trying to keep only the cells I previously kept after cleaning by mito and ribo thresholds in seurat saving into rds and exported in csv file...

But loom file come into the game and it's not fun.

Any help would be welcome. Thanks.

ldat <- ReadVelocity(file = "data.loom")
bm <- as.Seurat(x = ldat)

ids <- read.csv(file = 'daata_cellID_obs.csv', header = TRUE)

bm[["RNA"]] <- bm[["spliced"]]
dim(ids)
colnames(ids) <- "cells"
dim(bm)
bm <- subset(object = bm, cells = ids$cells)

It says "argument x is missing with no default value."
> [1] 11178 1 [1] 33694 12011 Error in UpdateSlots(object = x) :
> l'argument "x" est manquant, avec aucune valeur par défaut Calls:
> subset ... subset.Seurat -> UpdateSlots -> sapply -> lapply ->
> slotNames -> is Exécution arrêtée

Answer: Finally end up with the python way.

import numpy as np
import pandas as pd

from anndata import AnnData
from scanpy import read

import scvelo as scv

import loompy

# export HDF5_USE_FILE_LOCKING='FALSE'
# ssh -X

print("sample_obs")
sample_obs = pd.read_csv("/data/villemin/data/toulouse/scRNA-ALL/CellRanger/clusters/4006_rouge.4006_verte.CTL_Rouges.CTL_Vertes_cellID_obs.csv")
sample_obs.rename(columns={'Cells(seurat.Object)':'CellID'}, inplace=True )

sample_obs['CellID'] = sample_obs['CellID'].str[::-1]
sample_obs['CellID'] = sample_obs['CellID'].str.replace('_',':',1)
sample_obs['CellID'] = sample_obs['CellID'].str[::-1]
sample_obs['CellID'] = sample_obs['CellID'].str.replace('-1','x')

print(sample_obs.head())

filename = "/data/villemin/code/singleCell/python/4006_CTL_OSI.loom"
outputname = "/data/villemin/code/singleCell/python/4006_CTL_OSI.filtered.loom"

with loompy.connect(filename) as ds:
view = ds.view[:, np.isin(ds.ca.CellID, sample_obs['CellID'])]

loompy.create(outputname, view.layers, view.ra, view.ca)
Question: Hey everybody,
probably stupid question, but I do not understand what is wrong there, I guess I just do not know how to properly open a file in a function. Having my csv file in the folder, I want to make a list out of the it:

def make_OTU_list(csvfile):
with open(csvfile, "rb") as list_input:
csv = csv.reader(list_input)
list_OTU = list(csv)[0]
return list_OTU

make_OTU_list("OTU95_dominating.csv")

and I come back with that:

Traceback (most recent call last):
File "labbook.py", line 35, in
make_OTU_list("OTU95_dominating.csv")
File "labbook.py", line 31, in make_OTU_list
csv = csv.reader(list_input)
UnboundLocalError: local variable 'csv' referenced before assignment

Can somebody solve my problem?
Thanks a lot!

Answer: Look at the error message: "local variable 'csv' referenced before assignment"

1. Have you imported the csv library?
2. If yes, you should also use another variable name than csv to avoid confusion.
3. You probably want to write the list of the function into a variable.

( 4. Not really a bioinformatic question )

import csv

def make_OTU_list(csvfile):
with open(csvfile, "rb") as list_input:
csv_reader = csv.reader(list_input)
list_OTU = list(csv_reader)[0]
return list_OTU

your_list = make_OTU_list("OTU95_dominating.csv")
Question: I've a table of blast results with ~5hits/query protein:

```
Protein Class
ProtA 1
ProtA 1
ProtA 1
ProtA 0
ProtA 1
ProtB 1
ProtB 1
ProtB 0
ProtB 0
ProtB 1
```

I would like to convert this into a feature vector matrix like this:

```
Protein Class1 Class2 Class3 Class4 Class5
ProtA 1 1 1 0 1
ProtB 1 1 0 0 1
```

Can someone suggest me an efficient way to do this, since I've ~2300k hits in the file.

Answer: Using awk grouping

awk -F"\t" '{if(a[$1])a[$1]=a[$1]"\t"$2; else a[$1]=$2;}END{for (i in a)print i, a[i];}' OFS="\t" YOUR_FILE

outputs

ProtB 1 1 0 0 1
ProtA 1 1 1 0 1

I excluded the first/header line in the input. You can sort the rows if you want in any particular order.

Modified from: http://www.theunixschool.com/2012/06/awk-10-examples-to-group-data-in-csv-or.html
Question: Does anyone know how to suppress warnings from BioPython? Specifically, the Experimental warning from the use of SearchIO?

[This][1] script is perfectly functional, but I'd like to get rid of the annoying message about SearchIO - I'm using it for such a basic function that I'm not concerned that it's still experimental

[1]: https://github.com/jrjhealey/bioinfo-tools/blob/master/Genbank_slicer.py

I've tried:

>>> import warnings
>>> from Bio import BiopythonWarning
>>> warnings.simplefilter('ignore', BiopythonWarning)

as per the BioPython website, but that doesn't work in this case. Changing to BiopythonExperimentalWarning isn't correct either.

Answer: With some help from Peter Cock on twitter, it turns out that `warnings.catch_warnings():`, is the correct way to go, but it should be around the module import, not the function call:

# Import SearchIO and suppress experimental warning
from Bio import BiopythonExperimentalWarning
with warnings.catch_warnings():
warnings.simplefilter('ignore', BiopythonExperimentalWarning)
from Bio import SearchIO
Question: Hello,

I wanted to discuss frameshift mutations that difficult to interpret. I used both VCF file as well bam file to visualize them but I am unable to interpret what does bam or VCF telling me about that mutation.

e.g. in case of Bam file, IGV tells me that DEL 21 INS 6 G870(100%), T 0(0%)

Kindly can someone help me understanding such information? also if someone can explain the VCF information as well.

Thank you

Answer: Frameshift is a functional consequence of an insertion or a deletion, which creates an altered protein sequence with different amino acid content versus the reference protein (read [here][1] for more information). If you try to visualize a frameshift in IGV, you will only see the inserted or deleted nucleotides in the given position, and how many reads support that finding. Maybe uploading an screenshoot of what you need help with when visualizing the variant in IGV could help to have a more precise answer :).

You can find a lot of information regarding VCF format. You can always start with [wikipedia][2]. Basically each row in the VCF represents a variant position (unless you have a gVCF in which all the positions of the genome are listed) and the columns show different information, metrics, stats.. about the given variant.

[1]: https://www.genome.gov/genetics-glossary/Frameshift-Mutation
[2]: https://en.wikipedia.org/wiki/Variant_Call_Format