cannin/biostars_qa huggingface dataset:
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: Poisson distribution accounts for technical variance. NB distribution accounts for both technical and biological variance.

NB distribution is also a Poisson-gamma mixture distribution. Imagine you have a single biological sample (RNA extract) that you take aliquots out of to make technical replicates. These technical replicates will be Poisson distributed.

Now imagine you have multiple biological samples. You take multiple technical replicates out each biological replicate. You essentially now have multiple Poisson distributions for each biological replicate. The multiple Poisson distributions for each biological replicate can be described by a gamma distribution. Thus NB distribution (Poisson-gamma mixture) is used for RNA-seq.

You can also think of it as the lambda variable of the Poisson distribution is gamma distributed.
Question: Hi all,

I'm new to the community so apologise if this is a basic question or has been dealt with in another thread. I've looked through the fgsea documentation and am still having no luck.

I'm trying to run fgsea() on some nmf-decomposed gene expression data and keep getting an empty dataframe. This is confusing me because I conducted the analysis using MSigDB's GSEA software, and this seemed to work - I got rank-ordered gene lists and enrichment plots for each column of W_dash, so I'm sure it can be done in fgsea.

I ran nmf() on some gene expression data, then normalised the resulting W matrix. Now my task is to run GSEA on the normalised W matrix (W_dash), using the columns of W_dash as my ranked gene lists. W_dash looks like this:

W1 W2 W3 W4
AC000066_at 5.054266e-23 1.022037e-03 1.161830e-03 9.154324e-04
AC000099_at 2.620042e-03 1.369911e-03 1.495072e-02 5.842306e-03
AC000115_cds1_at 4.382477e-04 1.945546e-10 3.731523e-10 6.260872e-10
AC002045_xpt1_at 5.716331e-03 1.776428e-02 2.414694e-02 2.651467e-02
AC002073_cds1_at 3.328479e-03 4.652305e-04 4.209797e-03 4.125737e-03
AC002077_at 5.054266e-23 1.448357e-03 8.369223e-04 1.646869e-03
AC002086_at 1.109709e-03 1.874855e-03 3.724489e-04 2.646966e-03
AC002115_cds1_at 4.459879e-02 8.844502e-02 1.009927e-01 2.488297e-01
AC002115_cds3_at 2.441240e-02 1.509489e-02 2.094655e-02 1.398258e-02
AC002115_cds4_at 1.430019e-10 1.237761e-03 7.450935e-10 5.510702e-10
AC002115_rna2_at 8.419673e-03 1.567666e-02 7.554787e-03 3.770747e-02

I import a gene set from MSigDB:

>pathways = fgsea::gmtPathways("./c4.cgn.v7.4.symbols.gmt")

pathways looks like this:

>pathways[1]

$MORF_ATRX
[1] "ADCY3" "SEC31A" "BTD" "LTBP4" "UTRN" "FIG4" "AQP5" "TMEM11" "HNRNPL"
[10] "TAF5L" "INSIG2" "AGPS" "LSM1" "PIK3CB" "HTR4" "UBR2" "ZNF500" "RFC5"
[19] "ST13" "SSTR5" "PCF11" "CLP1" "MSX1" "CEP350" "COLQ" "HTR7" "MC2R"
[28] "NEK9" "PIAS2" "AMFR" "PDXDC1" "ZNF133" "MRE11" "MTX1" "CACNB2" "TMCC1"
[37] "PRELID3A" "JRK" "PPP1R3D" "NTPCR" "NKRF" "MLN" "KLHDC10" "OARD1" "FDXR"
[46] "PAXIP1" "WDR62" "RSRP1" "IPCEF1" "COX6A2" "POP4" "TCF12" "RERE" "LTN1"
[55] "ZNF292" "LEPROTL1" "KAZN" "TTLL5" "KRT33A" "GTSE1" "SEC62" "PARN" "MIA2"
[64] "DNAJC16" "SLC24A1" "SYNJ2" "SIK3" "MC5R" "FANCG" "GOLGA3" "MSL3" "ZBTB22"

Now I need to run fgsea() on the first column of W_dash as an ordered list:

>stats= sort(W_dash[,1], decreasing = T)
>stats[1:5]
M63379_at hum_alu_at S40719_s_at V00594_s_at Z70759_at
1.0000000 0.9656315 0.9164475 0.6934949 0.6814499

When I run fgsea() I get an empty dataframe:

>fgsea <- fgsea(pathways=pathways, stats=stats)
>fgsea
Empty data.table (0 rows and 8 cols): pathway,pval,padj,log2err,ES,NES...

Since the GSEA software produces the expected output, this leads me to believe I'm missing something really simple, but I can't work it out. Any help/guidance would be greatly appreciated!

Many thanks,

charlielonergan

Answer: The names of your ranked stats and the names of the entries of `pathways` must be the same. You have gene names in one but probe identifiers in the other. How should `fgsea` know the connection between both if the names are different?
Question: Hi everyone!

I am trying to find the best way to make 2 boxplot for a specific gene from data found in a row for a subset of columns within data frame x.

x dimensions are 634 by 128 columns

Each row is specific to a gene,

Column 1 has gene name, and I want to say look at gene in row#1

columns 2:48 data I want to include in one boxplot

columns 49:128 data I want to include in another boxplot

data frame looks something like this

```
gene accepted_hits_x1.bam accepted_hits_x1.bam etc....
1 AARS1 -6 0 etc....
```

I also want to be able to see each data point that makes up the boxplot plotted in the plot

I am having a problem:

I am running into the problem where my data (residual from mean ... meaning x value - mean) is a series of positive and negative values and it appears that with this plot it is excluding these negative values...

```r
data <- unlist(subset(datavr, gene =="IGF1R", select=2:128))

news <- data.frame(data=data, factor=c(rep(1,47), rep(2,80)))
news$data <- (log10(as.numeric(news$data)) + 1)

g <- ggplot(data=news, aes(x=as.factor(factor), y=data))

g + geom_boxplot() + geom_point(color="purple", size=3) + xlab("A38-41 A38-5 ") + ylab("log10(Residual from Mean)+1") + ggtitle("IGF1R inside region") + theme(plot.title = element_text(face="bold"))
```

The problem is that it keeps giving me error saying that:

```r
Removed 110 rows containing missing values (geom_point)
```

This could be that these values are negative so taking the log10(value)+1?

Answer: Since you added additional information, I'll just post this as an answer. Your best bet if you want to manipulate details about graphing in R is to use ggplot2 along with factors:

```r
library(ggplot2) # or install.packages("ggplot2"); library(ggplot2)

data <- unlist(subset(your_dataframe, gene == "geneName", select=2:128))
newData <- data.frame(data=data, factor=c(rep(1,47), rep(2,80))

g <- ggplot(data=newData, aes(x=as.factor(factor), y=data)
g + geom_boxplot() + geom_point(color="dark red", size=3) + xlab("x axis label") + ylab("y axis label") + ggtitle("My Plot Title") + theme(plot.title = element_text(face="bold"))
```

You can edit virtually everything you see with ggplot2; I only included the basics here. A google search for more will help with that.
Question: Hi,

I have a paired end DNA sequence alignment BAM file (got from using BWA -M) and would like to know if the points mentioned below is a good practice..

1) Since my aim is to perform variant calling, is it OK to keep unmapped, secondary,singleton and duplicated marked reads in the BAM file itself because by default they won't be considered for variant calling (bcftools mpileup and call)..OR..remove them to make the BAM file smaller to save space.

2) Since I am only interested in certain genes in the genome, I am planning to chop my BAM file to get only the reads aligned to my gene of interest (and 10000 bp upstream and downstream) and then perform variant calls on them. I can then perform variant calling on all my candidate genes of interest in parallel.

Can anyone let me know any caveat in the above mentioned approaches?

Answer: 1) I don't think the decrease in size would be a substantial space saver. I'm against all manual tampering with bam files. A good variant caller should ignore all non-proper reads.

2) You can do (at least in GATK) variant calling for specific intervals using the `-L` option. If you would start chopping bams you'll get in serious trouble with coordinates.
Question: Hi,

Currently getting the following error:

library(biomaRt)
ensembl <- useMart("ensembl")
Extra content at the end of the document
Error: 1: Extra content at the end of the document

Any ideas?

Cheers,

Answer: We're doing maintenance right now. Should be back in the next few minutes.
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: based on the first column only:

awk ' $1 !~ /@/ {print $1}' toy.sam | uniq -c

if you want to count only aligned reads you can also use awk to count only aligned reads.
Question: Hi, I am new to clusterProfiler and following the R documentation.

In the website instruction for enrichplot,
here is some demo for enrichGO:

library(clusterProfiler)
data(geneList, package="DOSE")
de <- names(geneList)[abs(geneList) > 2]
ego <- enrichGO(de, OrgDb = "org.Hs.eg.db", ont="BP", readable=TRUE)
library(enrichplot)
goplot(ego)

And I just use this but get no results.
Here is my code:

library(clusterProfiler)
library(org.Hs.eg.db)
library(DOSE)
data(geneList)
de <- names(geneList)[abs(geneList) > 2]
ego <- enrichGO(de, OrgDb = "org.Hs.eg.db", ont="BP", readable=TRUE)
head(summary(ego))

and I get:

> head(summary(ego))
[1] ID Description GeneRatio BgRatio pvalue p.adjust qvalue geneID Count
<0 rows> (or 0-length row.names)
Warning message:
In summary(ego) :
summary method to convert the object to data.frame is deprecated, please use as.data.frame instead.

So, what's wrong about my performance?
Could anyone give me a hint.

Answer: Hey!

Could you please show me with your contents of 'de', the input of your enrichGo command.

Because for my case, I will try to change my gene ID from 'ensembl_ID' to 'entrezid' to get more mapping results.

Since you got nothing mapping results, I guess there is something wrong in your input genes.

#Group_4 enriched Go terms

group4_top100_gene<- as.character((markergenes$Ensembl_ID)[9920:10019])

ID <- bitr(group4_top100_gene, fromType="ENSEMBL", toType="ENTREZID", OrgDb="org.Dr.eg.db")

head(ID)

ego <- enrichGO(gene = as.vector(ID$ENTREZID),
OrgDb=org.Dr.eg.db,
ont = "BP",
minGSSize = 1,
pvalueCutoff = 0.05,
qvalueCutoff = 1,
readable = TRUE)
head(summary(ego))

dim(ego)

result <- simplify(ego, cutoff=0.7, by="p.adjust", select_fun=min)

barplot(result, showCategory=10,title="GO terms enriched in group_4")

dotplot(result,showCategory=15,title="GO terms enriched in group_4")
dim(result)

Bests
Shizheng
Question: Dear All,
Am I right to use the following script to transform a fastq file (named test.fastq) to a fasta file? THANKS a lot!

#!/usr/bin/perl
use strict;
use Bio::SeqIO;
my $in=Bio::SeqIO->new(-file=>"test.fastq",-format=>'fastq');
my $out=Bio::SeqIO->new(-file=>">test.fasta",-format=>'fasta');
while(my $seq=$in->next_seq()){
$out->write_seq($seq);
}

Answer: This is the fastest bash one-liner to convert fastq to fasta, assuming 4 lines per FASTQ record, more on https://github.com/stephenturner/oneliners

sed -n '1~4s/^@/>/p;2~4p' test.fastq > test.fasta
Question: Hi Biostars,

I am looking for an open source software that has been widely used in the past and current with biologist, I am a software graduate student and my main focus is on scientific software testing. Greatly appreciate if you can name some of the open source commonly used software in your community?

Edit: So the main purpose of this query "Finding the most commonly used software" is to run some software tests on source code and see if we can help and improve open source software.

Answer: 1. BLAST

Have a look at the Biostars tags: https://www.biostars.org/t/ and try to extract only those tags which are software:

Then you get

1. R ahead of
2. Python
3. blast
4. samtools
5. Biopython

But R and Python are also programming languages (and have an 'unfair advantage') which you can do a lot of different things with, therefore blast is - not surprisingly - the tool Biostars' users have most questions about.
Question: I am doing a Preranked GSEA Analysis on the output of a differential expression analysis (genes ranked by sign(shrunken lfc) * -log10(pvalue) ).

I notice that when I feed the exact same files into GSEA, it can give me fairly different FDRs.

For example, when I ran the analysis the first time, it said it detected "32 gene sets are significantly enriched at FDR < 25%".

Ten minutes later, I ran the same files again, and it detected "4 gene sets are significantly enriched at FDR < 25%".

I know that there is some element of randomness in the GSEA algorithm because it sets a seed based on timepoint, but I'm surprised to see such a large difference in the number of pathways passing threshold from run to run.

Is this to be expected?

Thanks for your help.

**UPDATE** I also notice that the four gene sets called enriched in Round 2 are completely non-overlapping with the thirty-two gene sets enriched in Round 1. This is surprising and a little concerning to me.

Answer: It seems the problem is either that I needed to run more permutations or just that I had a very weak phenotype.

Pablo Tamayo at UCSD wrote the following in response to my question on the gsea-help listserv: " If you explicitly set the seed and get the same results then that is the expected behavior. The variability across different seeds (runs) is also expected to be higher on the FRDs than other statistics. The lack of robustness you are observing can be related with not running enough permutations and can be exacerbated because the results are rather weak to begin with. Have you tried running a much larger number of permutations?"

In response to his prompt, I reran the analysis with a higher number of permutations (initial # of permutations was 1000). The pathways called significantly enriched by FDR were still very variable, but the pathways called significantly enriched by nominal p-value much less so.

I decided eventually that this likely meant I was trying to chase down a very tiny-or-nonexistent phenotype, as I have since analyzed other databases with much more obvious and stable outcomes.
Question: Dear Users,

Recently I have aligned my Illumina reads to my reference genome. I have gone through sam file specifications but, it is quiet difficult to understand.

It would be great help if some one could explain the bam flags below.

A few lines of my aligned bam file:

```
HISEQ:115:C49JYACXX:4:1114:10901:17520 163 scaffold1 1 7 4M2I94M = 208 307 GCACAGAACTTTTTATTGAACAATTGGGCTACTAAATATGCTTACCAAACTCTATGGTAGTGTAGTCACAGTCTCTTAGTAGTAGTATCGATATATTGTA CC@FFFFFFHHHHJJJJIJJJIJJJIJJJGJIJIIJIJIIIJIJJGIJIIHIIJJJJFHIGHGBHCGIJIJIJJJGIJJCEHHHHHHFFFFDDDEFEEEE AS:i:-21 XS:i:0 XN:i:0 XM:i:2 XO:i:1 XG:i:2 NM:i:4 MD:Z:0A2G94 YS:i:-11 YT:Z:CP
```

Thank you

Answer: There is a neat app to explain the SAM flags: http://broadinstitute.github.io/picard/explain-flags.html
Question: I cannot decide which reference I should use to map my gut microbiome 16S data.
Greengenes seems to be the best choice, but last time it was updated in 2013 and metagenomics has moved far ahead since then. Won't I be missing a huge part of biodiversity, if I use it?

Also, this [article][1] keeps me worried. SILVA and Greengenes have a rather small overlap in taxonomy. It makes me want to make my own reference DB, but that'd be a serious standalone project.

What do you use to map gut 16S? Is there really nothing better than these two arguably accurate ref DBs?

[1]: https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-017-3501-4

Answer: Non-profit -> SILVA

For-profit -> RDP

Greengenes is dead, don't bother with it. Overall I think SILVA might be more accurate/up-to-date than RDP. However, SILVA has a costly license for non-academic settings..
Question: Hi all,

I am using DESeq2 for DE analysis and I want to normalize my read counts on the basis of two calibrator set (their counts) that were used. I am using estimateSizeFactors to give the read counts of calibrators as "controlGenes" so that I can run DESeq2 but when I execute my codes (which are as follows-)

calidata<-read.csv("/Users/pg2597/Desktop/combined/caliisocounts_xyz.csv", header=TRUE, row.names = 1)

mat_data <- data.matrix(calidata)
calibrators <- grepl('cali', rownames(mat_data))
cal <- calibrators
dds_norm <- estimateSizeFactors(mat_data,controlGenes=cal)

I get following error:

Error in (function (classes, fdef, mtable) :
unable to find an inherited method for function ‘estimateSizeFactors’ for signature ‘"matrix"’

Can someone help me with this?

Thanks in advance!

Answer: The `estimateSizeFactors()` function expects a DESeq data object, not a data matrix.

You should do this:

calidata<-read.csv("/Users/pg2597/Desktop/combined/caliisocounts_xyz.csv", header=TRUE, row.names = 1)

mat_data <- data.matrix(calidata)

calibrators <- grepl('cali', rownames(mat_data))

cal <- calibrators

dds <- DESeqDataSetFromMatrix(countData = countData, colData = colData, design = ~ condition)

dds_norm <- estimateSizeFactors(dds, controlGenes=cal)

For further information on this, including what the contents of *colData* should be and how to set the *design* parameter, see Quick start and Count matrix input.

Note that people usually use the wrapper function, `DESeq(dds)`, which will do the entire normalisation process for you, including calling `estimateSizeFactors()`. Your normalised counts could then be accessed with:

counts(dds, normalized=TRUE)

In this case, you are aiming to use control genes though.

Kevin
Question: Hello,when I download the .cel files from the [Broad CMAP][1] and open them, what I see is just digital gibberish.So I download the software called Gene Expression from the Aff company.

However,there are three algorithm named RMA,MAS5,PLIER WorkFlow to choose and the results are different when choosing different algorithm.

So which algorithm should I choose? And what difference are the three algorithm?

[1]: http://www.broadinstitute.org/cmap/cel_file_chunks.jsp

Answer: Run the RMA method.

See: http://bmbolstad.com/misc/ComputeRMAFAQ/ComputeRMAFAQ.html
Question: Hi all,

I've got some bulk RNA-seq data (mouse brain/hippocampus) and I'm looking for tools to extract celltype composition proportions from the data. I've been looking around and most of the tools focus on using associated single-cell data, which I do not have.

Do you know of any tool which provides celltype predictions (mouse brain) by only inputting bulk RNA-seq (e.g. using internal cell reference datasets)?

If not, is the best strategy to get some public single-cell data to use with my bulk RNA-seq (from what I've read, via MuSiC or CIBERSORTx for example)? In that case, could you recommend (if there exist) databases providing single-cell data with annotated celltypes (instead of raw, unannotated data) which I could use?

If not, what would be the best strategy to achieve this goal?

thanks for your help

Answer: >Do you know of any tool which provides celltype predictions (mouse brain) by only inputting bulk RNA-seq (e.g. using internal cell reference datasets)?

As soon as you search for tools published before 2019, the majority of deconvolution methods will probably be based on bulk RNA-seq alone.
See, for example, [this review from 2018](https://academic.oup.com/bioinformatics/article/34/11/1969/4813737).

>If not, is the best strategy to get some public single-cell data to use with my bulk RNA-seq (from what I've read, via MuSiC or CIBERSORTx for example)?

[Dorshizi et al, 2021](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8219045/) demonstrate the use of reference scRNA-seq data sets for brain-specific questions.

>In that case, could you recommend (if there exist) databases providing single-cell data with annotated celltypes (instead of raw, unannotated data) which I could use?

* The [Allen Brain Atlas](https://portal.brain-map.org/explore/transcriptome) is an incredible resource for transcriptomic (scRNA-seq) data from Mouse and Human data sets. I usually find a tad bit easier to first identify the paper with the corresponding data that I might be interested in and then track down the data on the Allen Brain Atlas servers, but I'm sure you can pick your way through.
* The [scRNAseq](http://www.bioconductor.org/packages/release/data/experiment/vignettes/scRNAseq/inst/doc/scRNAseq.html) package from bioconductor has a range of scRNA-seq data sets ready to load including cell labels
* the [SingleCell Portal[(https://singlecell.broadinstitute.org/single_cell) may be of help, too.

Apart from that, identifying publications with data sets of interest, e.g. via Pubmed or via GEO, should lead you to promising scRNAseq data set candidates; most of them should be provided with metadata (e.g. cell labels).
Question: Hi,
I just got a list of genes which were highly expressed in my cells.
I want to determine which genes are transcription factors.
Now I am searching the genes one by one but it needs a long time.
So is there a better and faster way?
Such as a database, where I input all the genes, then the transcription factors will be identified automatically.

Thanks!

Bing

Answer: Thank vchris_ngs very much. I have got the TF list in my gene list. Another question is that I want to do a network analysis. I have one target gene. I just want to know whether these TF are related with my target gene. Is there any good suggestions?
Question: Hey guys, I am doing de novo assembly. The data is downloaded from NCBI and they are DNA-seq short paired-end reads, 100bp x 2, about 120G in total. I use SOAPdenovo to assembly them. The k-mer value was set 29. After 10 hours, the server reported like this:

Ran out of memory while applying 13290854152bytes
There may be errors as follows:
1) Not enough memory.
2) The ARRAY may be overrode.
3) The wild pointers.

It was done on a server with 768G RAM and 50T internal storage. I don't know why it ran out of memory while applying 13G data, could anyone explain a little bit? There was no other tasks during the assembly.

Thank you.

Answer: Did you use the latest version? There was a bug giving the same error in previous versions, see [here][1], they seem to have fixed it in newer versions.

[1]: https://groups.google.com/forum/#!topic/bgi-soap/y0Q8g9WUFa0
Question: Hi,

I'm looking at some total RNA-Seq data and I'd like to get a quick look at the read distribution across "gene_type", according to Gencode annotations. Briefly, Gencode GTF files (at least for hg38) include annotations for the "gene_type", which identifies whether the region is protein coding, lincRNA, miRNA, etc. I'd like to know what percentage of reads originate from each type. RSeQC and RNA-SeQC seem to offer similar features (e.g., by exon, intron, UTR, etc.), but I don't think either provides by "gene_type".

Anyone know of existing software that will do this? I was planning to modify RSeQC read_distribution.py, but that works from BED file annotations (i.e., doesn't have "gene_type" annotations). RNA-SeQC doesn't provide source code. I could de-compile it and modify, but I don't see a specific license.

I appreciate your help.

Answer: You should be able to do this with [CGAT][1].

Briefly you want to prepare a bed file where the name column of the bed is taken from the gene_type attribute of the GTF. You can do this with

cgat gtf2bed --set-name=gene_biotype -I gencode.gtf.gz -S gencode_gene_type.bed.gz

You can then use `cgat bam_vs_bed` to count the number of reads that fall into each bed category.

cgat bam_vs_bed -a my_bam.bam -b genecode_gene_type.bed.gz > biotype_counts.tsv

[1]: https://www.cgat.org/downloads/public/cgat/documentation/cgat.html#cgat