cannin/biostars_qa huggingface dataset:
Question: I've got a list of SV calls, including breakpoints, and can easily enough winnow them down to those that are candidate gene fusions (intersect gene body or intron, same direction, strand, etc).

Now I'd like to know whether they're predicted to give in-frame or out-of-frame fusions. So far, I'm unable to find any annotation tool that can do this in a straightforward manner. Visualization would also be nice, but isn't necessary.

Any suggestions?

Answer: I had good results using [PRADA][1], (suggested by Roel Verhaak on Twitter). Specifically, the prada-frame command takes easy input (gene name, breakpoint location) that I already had and spits out a list of consequences in every transcript that matches.

[1]: http://bioinformatics.mdanderson.org/main/PRADA:Overview
Question: Hi

I've been working on the gene set enrichment analysis by ClusterProfiler, and I'm able to generate the enriched terms table by the `enricher` function:

CRE_MapMan_GSEA <-enricher(DEG_cht7_6_up_LocusID$Gene_ID,
TERM2GENE=TERM2GENE,
TERM2NAME = TERM2NAME,
pvalueCutoff = 0.05)
Result:

ID	Description	pvalue

28.1	DNA.synthesis/chromatin structure	6.96E-19
31.6.1.1	cell.motility.eukaryotes.basal bodies	4.77E-08
28.2	cell.motility.eukaryotes.basal bodies	1.63E-06
24	DNA.repair	4.38E-06
31.2	Biodegradation of Xenobiotics	1.32E-05
27.3.44	cell.division	1.90E-03
31.1	RNA.regulation of transcription. 2.24E-03
Chromatin Remodeling Factors

However, there're only 7 terms displayed. Does that mean the rest of them are not significant enough? Is there any way I can see the p-values of more terms? I've tried to play with the `pvalueCutoff` parameter, but the number of the displayed term is still the same.

Thanks!

Answer: Try to use different `qvalueCutoff` value as cut off.
Question: I am trying to convert these gene symbols( that's what i think they are) to ensembl ID's using biomaRt package. But for some reason it seems to fail. Would love to have some feedback!

library(biomaRt)

ensembl <- useMart("ensembl", dataset="mmusculus_gene_ensembl")

genes <- c("Acad11","Acad8", "Acad9","Acadm","Acadsb", "Acan" ,"Acap1","Acap2" ,"Acat1","Acbd3", "Acbd4")

foo <- getBM(attributes=c('ensembl_gene_id',
'external_gene_name'),
filters = 'hgnc_symbol',
mart = ensembl,
values = genes)

Thanks in advance.

Answer: Almost correct but it must be `mgi_symbol` instead of `hgnc_symbol` since the HGNCs represent human genes. MGI is the equivalent for mouse.
By the way it is **Ensembl** not Ensembl**e**, just saying if you ever write this up as a paper be sure to correctly spell the data sources ;-)

library(biomaRt)

ensembl <- useMart("ensembl", dataset="mmusculus_gene_ensembl")

genes <- c("Acad11","Acad8", "Acad9","Acadm","Acadsb", "Acan" ,"Acap1","Acap2" ,"Acat1","Acbd3", "Acbd4")

foo <- getBM(attributes=c('ensembl_gene_id',
'external_gene_name'),
filters = 'mgi_symbol',
mart = ensembl,
values = genes)

> foo
ensembl_gene_id external_gene_name
1 ENSMUSG00000090150 Acad11
2 ENSMUSG00000031969 Acad8
3 ENSMUSG00000027710 Acad9
4 ENSMUSG00000062908 Acadm
5 ENSMUSG00000030861 Acadsb
6 ENSMUSG00000030607 Acan
7 ENSMUSG00000001588 Acap1
8 ENSMUSG00000049076 Acap2
9 ENSMUSG00000032047 Acat1
10 ENSMUSG00000026499 Acbd3
11 ENSMUSG00000056938 Acbd4
Question: Given a sam file which was produced by using bwa to align a fastq of reads to a reference genome, how can I enumerate the number of reads which were assigned to each chromosome of the reference genome?

Answer: Your question refers to SAM format, but in the event that you may also have a sorted/indexed BAM, you can simply do the following to get this very quickly.

samtools idxstats aln.bam | cut -f 1,3
Question: Hi,
I want to do psiblast in command line. Before, performing psiblast, I tried to make a blast database by using makeblastdb command. But it shows the following error:

BLAST Database creation error: Error: Duplicate seq_ids are found:
REF|WP_003261842.1
My fasta file is very large. So, manual deletion of the duplicate sequence is nearly impossible. How can I get rid from the duplication?

Cheers

Answer: Use [seqkit][1]:

$ seqkit rmdup input.fa > output.fa

[1]: https://bioinf.shenwei.me/seqkit/usage/#rmdup
Question: Hi guys,

I am trying to download RNA seq using fastq dump. The reads are paired end reads so I am also splitting them. When I try to download the 24 files, the sra cache for each file is also being downloaded in the home directory where I only 13 gb of space. The first couple of files download in the directory where I have 10TB of space but their .sra.cache files download in my home directory and it fills up space quickly and once it reaches the 13 gb, the rest of the files download but I get the error "No space left on device". I am wondering even though I see the files I downloaded, is something wrong with them or are they fine and why do I keep getting the error? Appreciate any help. Thank you.

Answer: Did you run the `vdb-config -i` to designate a location where SRA files are downloaded? Additional help at: https://github.com/ncbi/sra-tools/wiki/03.-Quick-Toolkit-Configuration
Question: I have this software I had to use for my PhD that was nowhere near as speedy or memory-efficient as I needed with the mountains of data we have. Therefore I rewrote it from scratch, outsourcing a lot of it to C and made it parallel (it was trivially parallelizable up to the number of chromosomes in the organism). Since it is so memory efficient it also allows for computations that were not previously possible.

1) Now, is this thing publishable? I have made no improvements to the results of the algorithm per se it is just improved from a user / comp.sci. perspective.

2) And what are the ethics of this? If it were published, I'd imagine I'd be getting the cites, not the original authors that actually invented the technique/algorithm. I'd feel bad about this part. My changes are mostly in the implementation, great as they may be.

Thoughts on these issues welcome. Same goes for any related issues I did not consider.

Ninja edit: 3) what are the pros and cons of contacting the original authors? I am afraid they'd just claim that I did not do anything original, tell me to publish my code on github and redirect the cites to their original paper.

Answer: 1. Sounds like you a software paper or application note or whatever you journal of choice calls short papers describing the implementation of some method.
2. Conscientious researchers will cite the "method of [those other people], as implimented in [software] (Brakken Stovner, 2015)". Since you will write your software's documentation you can explicitly ask users to do this. If your software makes the method more applicable you both get more cites.
3. I would contact them -- for one thing, they are the people most likely to be asked to review the software paper anyway. If you contact them you get to explain your motivation and how your software helps them, and you might get feedback from them that helps you software.
Question: Given a sam file which was produced by using bwa to align a fastq of reads to a reference genome, how can I enumerate the number of reads which were assigned to each chromosome of the reference genome?

Answer: The solution based on samtools idxstats aln.bam should be used with caution. This is because AFAIK the numbers reported by samtools idxstats (& flagstat) represent the number of alignments of reads that are mapped to chromosomes, not the (non-redundant) number of reads, as stated in the documentation. I stumbled across this by observing that the total numbers reported by samtools idxstats (& flagstat) exceeded the total number of reads in my input fastq files. See further discussion on this issue in the comments on this post: http://left.subtree.org/2012/04/13/counting-the-number-of-reads-in-a-bam-file/
Question: Hello everyone!

I tried to make a basic scatter plot with R using gene expression data.

```r
#import data:

oldmice <- read.table("oldmice.txt", header = TRUE)
youngmice <- read.table("youngmice.txt", header = TRUE)
```

Imported data contains: format is the same for both imported data but MGE has different values.

gene	MGE
Sox17	-6.74193774617653
Mrpl15	-0.212567471203473
Lypla1	-0.711251006455475

and so on..

Made basic volcano plot using: youngmice$MGE vs oldmice$MGE

```r
plot(oldmice$MGE, youngmice$MGE, main="old vs young mice!!",
xlab="oldmice$MGE ", ylab="youngmice$MGE ", pch=19)
```

My question is how to color "genes" which is in `multiple_gene_lists` into `oldmice$MGE`, `youngmice$MGE`? (which should label the only `multiple_gene_list` which are in `multiple_gene_lists` into `oldmice$MGE`, `youngmice$MGE`)

Here is my multiple_gene_list

```r
multiple_gene_list <- read.table("multiple_gene_list.txt", header = TRUE)
multiple_gene_list <- as.vector(multiple_gene_list )
```

`multiple_gene_list` contains:

gene
Six6
Arl2
Tmem74B
Rab9B
Rasgef1B
Ccne1
Apln
Spag7
C17Orf59
Krtap4-4

And my goal is to only label `multiple_gene_list` in `oldmice$MGE`, `youngmice$MGE`. I also tried the following code but failed!

```r
with(subset(ASC_oldmice_exprs, ASC_oldmice_exprs$gene %in% multiple_gene_list$gene), points(ASC_youngmice_exprs$MGE, pch=20, col="red"))
```

Thank you!

Answer: Hey, here is a simple way to do it, and to also have a different shape for the genes of interest:

oldmice <- data.frame(gene = c("BRCA1","BRCA2","BRCC3","TP53","ATM"),
MGE = c(-6,-5,-8,-7,-6))
youngmice <- data.frame(gene = c("BRCA1","BRCA2","BRCC3","TP53","ATM"),
MGE = c(-6.5,-7,-4,-7,-8))
multiple_gene_list <- c("TP53","ATM")

1, Create a default colour vector:
----
mycol <- rep("grey50", length(youngmice$gene))
mycol
[1] "grey50" "grey50" "grey50" "grey50" "grey50"

2, Set a colour for the genes of interest:
----
mycol[which(youngmice$gene %in% multiple_gene_list)] <- "firebrick1"
mycol
[1] "grey50" "grey50" "grey50" "firebrick1" "firebrick1"

3, Same as 1 & 2 but for shape (pch):
----
mypch <- rep(19, length(youngmice$gene))
mypch[which(youngmice$gene %in% multiple_gene_list)] <- 18
mypch
[1] 19 19 19 18 18

4, Plot data:
----
plot(oldmice$MGE, youngmice$MGE, main = "old vs young mice!!",
xlab = "oldmice$MGE ", ylab = "youngmice$MGE ", pch = mypch, col = mycol, cex = 5.0,
xlim = c(-10, -2), ylim = c(-10, -2))

text(
subset(oldmice, gene %in% multiple_gene_list)$MGE,
subset(youngmice, gene %in% multiple_gene_list)$MGE - 0.45,
lab = subset(youngmice, gene %in% multiple_gene_list)$gene,
cex = 1.0)

You can also set the labels this way:

mylab <- as.character(youngmice$gene)
mylab[-which(youngmice$gene %in% multiple_gene_list)] <- ""
mylab
[1] "" "" "" "TP53" "ATM"

plot(oldmice$MGE, youngmice$MGE, main = "old vs young mice!!",
xlab = "oldmice$MGE ", ylab = "youngmice$MGE ", pch = mypch, col = mycol, cex = 5.0,
xlim = c(-10, -2), ylim = c(-10, -2))

text(oldmice$MGE, youngmice$MGE - 0.45,
lab = mylab,
cex = 1.0)
Question: In all RNA-seq analysis applications they talk about the dispersion of a gene. As far as I understood, it is not a variance of the normalized counts for a given gene. It is somehow much more complicated.

DESeq defines the dispersion as follows:

But what would a dispersion of 0.19 or a dispersion of 0.80 tell me? Can I still interpret it as a variance of a gene?

Answer: https://support.bioconductor.org/p/75260/
Question: Hi everyone,

I would like to download some ENCODE data. I want to download which genes present in their promoter the binding of particular Transcription Factor. Starting from a trasncription factor I would like an easy way to retrieve a gene list (not sequences, but more processed data) where they find that Transcription Factor. The best, will be a way to access this data through R.

Thank you for your help

Answer: If it is human, try the [ReMap][1] database and intersect your regions with it.

[1]: http://tagc.univ-mrs.fr/remap/
Question: Hello everybody,

I am trying to analyse CNV data from TCGA to get a measure of overall CNV per patient.

When I download the Level 3 files taken from the SNP6 array, there is a column in the file called Segment_Mean. (Example at bottom.)

What do the numbers in this column represent?

I think they might be log ratios, but the link below makes me wonder if they are direct estimates of copy number. (In which case it is puzzling that they aren't whole numbers.)

http://www.broadinstitute.org/cancer/software/genepattern/modules/snp6copynumberpipeline

Thanks for any help,

Stephanie

```
Sample Chromosome Start End Num_Probes Segment_Mean
BONZE_p_TCGAb56_SNP_1N_GenomeWideSNP_6_E04_666936 1 151040529 153927851 1558 0.2031
BONZE_p_TCGAb56_SNP_1N_GenomeWideSNP_6_E04_666936 1 153928595 153929981 2 -2.0772
BONZE_p_TCGAb56_SNP_1N_GenomeWideSNP_6_E04_666936 1 153933585 164456865 7473 0.1883
```

Answer: Those are the log2 ratio of the tumor intensity to the normal intensity. To convert to an absolute cn, use: `(2^seg_mean)*2`
Question: I am running GWAS analyses with gemma that seem to work but when I open the centered relatedness matrix (.cXX), I do get a symmetric matrix but the main diagonal is not made of 1s ... why?

The [manual][1] p.11 part 3.3.1 shows an example of matrix with a main diagonal which is not made of 1s.
I also tried to dig into the [gemma code][2] on Github but I did not find the piece of code used to calculate the relatedness matrix.

---
####Commands and output

Plink and gemma commands to generate the matrix:

# Generate plink file from my VCF file (67 samples)
vcftools --gzvcf $vcf_file --plink --out $prefix_vcf
# Make binary plink files
plink --file $prefix_vcf --make-bed --out $prefix_vcf
# Generate centered relatedness matrix with gemma
gemma -bfile $prefix_vcf -gk 1 -o $prefix_vcf

Here the output log of gemma:

##
## GEMMA Version = 0.94
##
## Command Line Input = -bfile subset_67_accessions_wo_singletons_only_alt_allele -gk 1 -o subset_67_accessions_wo_singletons_only_alt_allele
##
## Summary Statistics:
## number of total individuals = 67
## number of analyzed individuals = 67
## number of covariates = 1
## number of phenotypes = 1
## number of total SNPs = 927314
## number of analyzed SNPs = 39678
##
## Computation Time:
## total computation time = 0.123833 min
## computation time break down:
## time on calculating relatedness matrix = 0.00733333 min
##

I get as expected a square matrix of order 67 (as my number of samples in the VCF file).
Here are the first 5 rows and columns of the matrix:

0.5265889259 0.0918289748 0.01315201643 -0.02417258371 -0.00300349980
0.0918289748 0.4426259465 0.01751500729 -0.0196104412 -0.04093555612
0.01315201643 0.01751500729 0.494663962 0.01773204554 0.03579284737
-0.02417258371 -0.0196104412 0.01773204554 0.5622735808 0.05660374681
-0.003003499807 -0.04093555612 0.03579284737 0.05660374681 0.6064086438

There is symmetry but not a main diagonal of 1s while same individuals should be identical to themselves, shouldn't they?

*NB*: I posted a similar post on the [gemma-discussion][3] google group but still did not get an answer after a month.

[1]: http://www.xzlab.org/software/GEMMAmanual.pdf
[2]: https://github.com/genetics-statistics/GEMMA/tree/master/src
[3]: https://groups.google.com/forum/#!topic/gemma-discussion/34FyII8oXUM

Answer: I've got the answer from the creator of GEMMA Xiang Zhou:

> The centered relatedness matrix does not scale each column of the genotype matrix, so the main diagonal of the resulting relatedness matrix will not be close to 1s (usually around 0.3 depending on maf etc.). If you use the scaled relatedness matrix (.sXX), then the main diagonal should be made of values that are very close to 1s.
Question: I’m not sure how familiar folks are with the Bioconductor ">Biobase" datasets but one of the free datasets is called ALL (Acute lymphoblastic leukemia). My goal is to calculate and visualize some basic properties of gene expression levels from the ALL dataset. I’m getting stuck on how to write a line of code that's specifically for patients, and another just for genes. I understand that rows=genes and columns=patients, but I do not have the syntax down to specify patients and genes. For example, I have tried:

> hist(exprs(ALL))
And

> median(exprs(ALL))

[1] 5.469092

I’m not happy with the results and I know I’m doing something wrong, for example…

> hist(exprs(ALL))

…will not calculate the average gene expression levels for each gene.
and…

> median(exprs(ALL))

…will not find the average gene expression per gene. **Any guidance is super appreciated!**

Answer: How about `apply(exprs(ALL), 1, median)`.
`apply` calls a function (here median), on each row (1) or column (2) of x, where x is `exprs(ALL)`.
Question: I am new to plink, and I am following [this tutorial](http://zzz.bwh.harvard.edu/plink/tutorial.shtml#t2):

But I cannot even open the file.

I am using a PLINK 1.9 on Windows, 64 bit. I have put the hapmap1.ped and hapmap1.map in the folder which I store PLINK.

I used the below command to add the folder to the path

path = %PATH%;C:\Users\jh3903\Desktop\plink_win64

And then when I run the below command, it return my PLINK version is `v1.9b04.5 64-bit (25Jul2017)` and with some other information

plink

I then run the below command

plink --file hapmap1

This is the error I got says

> 16306 MB RAM detected; reserving 8153 MB for main workspace.
> Error: Failed to open hapmap1.map

Does anyone have a clue about this?

[1]: http://zzz.bwh.harvard.edu/plink/tutorial.shtml#t2

Answer: Does the file `hapmap1.map` exist in the directory you're running the `plink` command in? It should be in your working directory, not in the directory that has the plink executable.
Question: Hi all,

I want to do read normalization using trinity, actually with looking the written command in the trinity's site, I'm confused about the allocated RAM to Jellyfish and coverage. It says:

```
TRINITY_RNASEQ_ROOT/util/insilico_read_normalization.pl \
--seqType fq \
--JM 100G \
--max_cov 30 \
--left left.fq \
--right right.fq \
--pairs_together \
--PARALLEL_STATS \
--CPU 10
```

Could anybody help me to know what `--JM 100G --max_cov 30` means. From what I understand, 100GB RAM is allocated to Jellyfish that seems is high for computer with 64GB RAM. please let me know the exact command for read normalization in my computer (64GB RAM). Thank a lot in advance

Answer: If you want to normalize reads and assemble the transcript with normalized reads (default 50x coverage):

TRINITY_RNASEQ_ROOT/Trinity --normalize_reads --seqType fq --JM 60G --left left.fq --right right.fq --CPU 6

If you want to change the parameter, say `--normalize_max_read_cov` to `30`, add the following to the above command: `--normalize_max_read_cov 30`
Question: Hi everyone!

I am a beginner using GATK, so bear with me please.

I am trying to do my variant calling with GATK4 new Mutect2 (not MuTect2), using as --germline-resource the af-only-gnomad.hg38.vcf.gz from the GATK's bundle. The command is pretty much the same as seen in the tutorial: https://software.broadinstitute.org/gatk/documentation/article?id=11136, only updated to allow it to be used in a WES.

My problem is that after successfully completing the variant calling, none of the variants have ID information, although when I compared it to a MuTect2 result several did. Basically, all of them have a "." in that field.

I checked the uncompressed gnomad.hg38.vcf.gz manually and the ID info as rsXXXXXX are there.

Any ideas on why the variant calling is not collecting the ID info from the germline resource to populate the ID field?
Thank you!

Daiana

Answer: Just in case someone encounters the same issue, it has been answered in: https://gatkforums.broadinstitute.org/gatk/discussion/comment/57540#Comment_57540

There is no way of populating the ID field with GATK4 Mutect2 as of today (03/29/19).

Best,

Daiana
Question: Dear Friends,

I have a matrix like this:

Gene BRCA THYM TGHJ
ACC 23 21 7
XTG 12 13 9
CFG 45 4 8

The numbers are the snp count from vcf file.

I want to plot this in the form of a matrix with numbers colored based on their count; for example the highest number is colored in "Red" and then gradually the intensity of the color decreases with the decrease in number, so in this 45 is colored "Red" and 4 is colored with a very light color. Please let me know if am clear.

I am looking to plot this matrix using ggplot2, but other ways in R are also very welcome. A matrix like this:

![enter image description here][1]

[1]: https://s8.postimg.cc/65y00ijs5/Screen_Shot_2018-08-02_at_1.05.39_AM.png

Answer: library(tidyr)
df2=gather(df1, Symbol,Count, -Gene)
library(ggplot2)
ggplot(df2, aes(Gene,Symbol, size = Count, label = Count)) +
geom_text(size = df2$Count, aes(colour=df2$Count))+
scale_colour_distiller(palette="RdYlGn")+
theme_bw()
Question: Dear All,

I am a Biologist trying to understand the statistics of RNA-Seq data.

Given that RNA-Seq follows NB distribution with Biological replicates, as NBD accounts for overdispersion in the data, I am not sure how to ascertain it to my data.

Although I understanstood these distributions through standard books I am unable to comprehend and relate it to RNA-Seq.

Differential expression is my aim.

I have simulated data, to test and understand some open source tools like edgeR, DESeq, Cufflinks etc.

I have real data set too.

I have two conditions with four replicates each.

If I have to know whether my data fits NBD or Poisson distribution, I have to check this across replicates of each gene of each condition??

If the above point is right, how do I do it?

Should I do some goodness of fit test like Chi-sqare test or just the mean variance relationship is enough?

Thanks in advance for your valuable inputs.

Answer: If you have biological replicates, then they're pretty much guaranteed to fit a negative-binomial distribution better than a Poisson distribution (otherwise, there's no biological variance). If you wanted to check, graph variance vs mean. If the values don't cluster on the dispersion==mean line, then it's not Poisson.
Question: I have 50bp PE RNA-seq, and 4 of my ~200 samples end in error after running Cufflinks (v.2.2.1). Here is my pipeline prior to running Cufflinks:

1. STAR 2-pass (GRCm38)
2. Added read group information with Picard
3. Marked (but did not remove) duplicates with Picard
4. SplitNCigarReads with GATK
5. Realigned around indels with GATK
6. Performed base recalibration with GATK

I've been looking at one particular sample with the error and have narrowed it down to a couple reads (see below):

PC168224:108:C90EBANXX:5:1303:19181:37569	179	6	85461128	60	43M345H	=	85461508	338	CCCTCAGATCATCATCCGAGCTTTCCGCACAGCCACCCAATTG	C46FE@3CDECDBDD>/C6FE?FD@:EA?5CDBB>33<::;

Answer: this is a bad cigar string:

380H4I4S

there is not "Match/M" in this cigar string: 380 hard clip, 4 insertions, 4 soft clip. Just 3 things that are not associated to a reference;
Question: I am interested in Bioconductor and I want to dive into their stats pages.

- [Index of the most used Bioconductor software packages][1]
- [The most used Bioconductor annotation packages][2]
- [The most used Bioconductor experiment packages][3]

Is there a way to get a downloadable file with the stats or do I need to do it one by one?

Is there any other parsable stats page? [This][4] is the only type of stats page

PS: Related to the last answer of [this post][5]

[1]: http://bioconductor.org/packages/stats/index.html
[2]: http://bioconductor.org/packages/stats/data-annotation.html
[3]: http://bioconductor.org/packages/stats/data-experiment.html
[4]: http://bioconductor.org/packages/stats/bioc/Biobase.html
[5]: https://www.biostars.org/p/60334/#60367

Answer: As this seems to be popular since then there is a [package][1] BiocPkgTools, currently only in github to download and analyse the data

[1]: https://seandavi.github.io/BiocPkgTools/articles/BiocPkgTools.html#download-statistics