cannin/biostars_qa huggingface dataset:
Question: I was seeking the pathway analysis of a gene list and just saw two contrast things which made me too confused.

In https://support.bioconductor.org/p/54827/ the guys discuss enrichment analysis and "Paul Shannon" has suggested the below code to calculate p-values.

```r
library(KEGGREST)
library(org.Hs.eg.db)

# created named list, eg: path:map00010: "Glycolysis / Gluconeogenesis"

pathways.list <- keggList("pathway", "hsa")

# make them into KEGG-style human pathway identifiers
pathway.codes <- sub("path:", "", names(pathways.list))

# for demonstration, just use the first ten pathways
# not all pathways exist for human, so TODO: tryCatch the
# keggGet to be robust against those failures

# subsetting by c(TRUE, FALSE) -- which repeats
# as many times as needed, sorts through some
# unexpected packaging of geneIDs in the GENE element
# of each pw[[n]]

# genes.by.pathway <- sapply(pathway.codes,
# function(pwid){
# pw <- keggGet(pwid)
# pw[[1]]$GENE[c(TRUE, FALSE)]
# })
load(paste(Data_path, "KEGG_Gene_Pathways.RData"))

all.geneIDs <- keys(org.Hs.eg.db)

# chose one of these for demonstration. the first (a whole genome random
# set of 100 genes) has very little enrichment, the second, a random set
# from the pathways themsevles, has very good enrichment

genes.of.interest <- c("23118", "23119", "23304", "25998", "26001", "51043",
"55632", "55643", "55743", "55870", "7314", "56254",
"7316", "144193","784", "8837", "1111", "84706",
"200931","169522","5707", "5091", "5901", "55532",
"9777")# the hypergeometric distribution is traditionally explained in terms of

# drawing a sample of balls from an urn containing black and white balls.
# to keep the arguments straight (in my mind at least), I use these terms
# here also

pVals.by.pathway <- sapply(names(genes.by.pathway),
function(pathway) {
pathway.genes <- genes.by.pathway[[pathway]]
white.balls.drawn <- length(intersect(genes.of.interest, pathway.genes))
white.balls.in.urn <- length(pathway.genes)
total.balls.in.urn <- length(all.geneIDs)
black.balls.in.urn <- total.balls.in.urn - white.balls.in.urn
total.balls.drawn.from.urn <- length(genes.of.interest)
dhyper(white.balls.drawn,
white.balls.in.urn,
black.balls.in.urn,
total.balls.drawn.from.urn)
})

print(pVals.by.pathway)
```

As you know, dhyper returns

On the other side, http://www.tongji.edu.cn/~qiliu/help/help_3_ORA.html provides a complete definition of pathway over-representation analysis, and says the p-value is calculated as

Actually, my problem is not about the summation over 'x', but it is about the subtraction of the formula from '1'.

Which of them is true? Does anyone know about it?

Answer: It's shown in the notation. You have to read the whole paragraph to understand the difference.

The first formula is the probability that X=k. That's not a hypothesis test p-value at all. If you want to know how uncommon such an event is, you would sum the probabilities of the more unlikely events, and take the difference from 100%. That's what's happening in the second formula.

In the second formula it's adding up the probabilities of some events, and showing the difference from 100%, or the rareness instead of the likeness.
Question: Dear All,
I have two networks. Both the networks have common nodes as well as unique nodes. Is there any way I can combine both and analyze??

Answer: Using **R** igraph union:

library(igraph)

net1 <- graph_from_literal(A-B:C:D, A-C, B-C)
net2 <- graph_from_literal(A-X:Y:Z, Z-Y)

par(mfrow = c(1, 3))
set.seed(1)
plot(net1, vertex.size = 50, main = "net1")
plot(net2, vertex.size = 50, main = "net2")
plot(union(net1, net2), , vertex.size = 50, main = "union of net1 and net2")
Question: Is there any way to convert a cytoband location (e.g., 17p13.3) to its genomic coordinates programmatically (i.e., without having to do a manual conversion)? A web API such as BioMart would be ideal, but I haven't found a way to do it this way. Entrez Eutils would also be okay, but I likewise haven't been able to figure out how to do it this way.

Answer: $ CHR="chr17"
$ LOCATION="p13.3"
$ wget -qO- "http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/cytoBand.txt.gz" | gunzip -c | grep -F "^${CHR}" | grep -F ${LOCATION} > answer.bed
Question: Deal all,
I exported human-pig ortholog relationship information from ensembl biomart. For gene SIRPA, it showed empty information of pig gene/transcrip/protein id nor any ortholog alignment information. But for most other genes, there are such information. So I thought pig doesn't have SIRPA. However, I found pig's SIRPA information when searching the gene in NCBI-gene. Why? I wondered maybe it's because they shared the same gene name but the properties of pig's SIRPA has changed quite a lot from human's SIRPA. But the gene description showed that they are quite the same. Could anyone explain this phenomenon a bit? Does ensembl biomart provide ortholog information incomplete? Thank you.

Best,
Yingzi

Answer: Orthologues in Ensembl are assigned using a [protein tree pipeline](http://www.ensembl.org/info/genome/compara/homology_method.html) which clusters proteins sequences using HMMs, then builds multiple alignments and reconciles against the species trees. The relationships between genes in these trees are used to [assign homologues](http://www.ensembl.org/info/genome/compara/homology_types.html).

These are recomputed with each Ensembl release and the pipeline was recently changed from using Blast for the clustering to HMMs, which has altered a lot of trees. The [previous version of Ensembl](http://jul2018.archive.ensembl.org/Homo_sapiens/Gene/Compara_Ortholog?anc=12844285;db=core;g=ENSG00000198053;g1=ENSSSCG00000007181;r=20:1894167-1940592;t=ENST00000356025) had human/pig orthologues with [many-to-many relationships](http://jul2018.archive.ensembl.org/Homo_sapiens/Gene/Compara_Tree?anc=12844285;db=core;g=ENSG00000198053;g1=ENSSSCG00000007181;r=20:1894167-1940592;t=ENST00000356025), which has changed. I think this needs to be re-examined and the orthologues restored. I recommend using the archive version of Ensembl for now.
Question: Hello everyone,

I'm pretty new to the field, so forgive me if this is a dumb question.
So, I am working on mapping out single reads to the reference mouse genome, with given 260 fastq files.
What I'd like to do is make all 260 files to be good quality reads. Right now it's a total mess, as in some have bad 'per base sequence quality' and 'per base sequence content', while some reads have bad 'per base GC content', and so on.

My original thought was run a program that fetches what passed and what failed from summary.txt from each fastqc directory. Then I'd have list of files that failed 'per base sequence content' or failed 'per base GC content' and such. However, because each reads are different within the list, I would still have to open html file and look individually to improve the quality. ( ex. I might have to trim 10 bp while other files need to trim more than that to improve the quality. )

Is there a way to somehow program it so that I can make all the files to be good quality reads at once, rather than checking each fastqc.html files myself and make different adjustments to each reads to improve the quality?

Answer: Try [multiqc][1], it collects all the fastqc reports into one overall report - extremely handy if you have lots of files.

[1]: http://multiqc.info
Question: I am working on wheat genome, I want to do analyze comparative genome analysis of 3 varieties of wheat .I have sequence files from illumina 1.9 in fastq format. I checked the quality of reads by fastqc tool.GC content are not in normal range.(47- 49). What is the normal value of % GC for RNA seq reads ? The other question is that kmers are also not in correct range.How can I correct it.
For trimming,adaptor seq file is required,but i don't have this file..Is that possible to remove these two error?If yes then how can I do?Can I skip the trimming step and go to next step of mapping?

In this file all parameter values are correct except kmer s and GC content.Is there any need to trim it?If yes then how can I do?

file:///home/comsats-ra/fatimamphilldata/G1_cleaned_R1_fastqc.html#M11

Answer: You can use [fastp](https://github.com/OpenGene/fastp) to trim adapters for Illumina sequencing data, without the need of knowing the adapter sequences.

Just download fastp and run:

fastp -i in.fq -o out.fq

And then everything is done, the adapters are trimmed in out.fq

For paired end data, the command is like:

fastp -i in1.fq -o out1.fq -I in2.fq -O out2.fq

Gzip is supported for both input and output.
Question: Hi all,

I have a VCF file that was generated using HG18 and I want to convert it to VCF but that has HG19 data.

I am aware that there is the lift-over tool from UCSC but how do I convert a VCF file from HG18 to a VCF file to HG19, since it seems the lift-over tools produces a BED file instead.

Thanks in advance.

Answer: You may use GATK "LiftoverVariants". It allows you to use VCF file as an input and outputs a VCF file corresponding to the new assembly. Link: https://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_gatk_tools_walkers_variantutils_LiftoverVariants.php
Question: Hello clever community!

I need your advice.
I am working on a *de novo* plant genome assembly of ~400 Mb.
I have Chromium 10x data, which was assembled with supernova. I also have Illumina paired end reads. Now I have additional data of PacBio reads, 120x roughly. The genome is diploid and I am thinking about using Falcon.

What do you think should be the best strategy:

1. Assembling PacBio reads and then using a tool to integrate the two assemblies? Is there anything like this? Which tool would you use?

2. Using a tool that can assemble the genome from both the chromium and the PacBio reads? Is there anything like it?

3. Assembling the PacBio reads and using chromium 10x and the illumina for polishing? If I assemble with Falcon, what tool should I use for polishing?

4? Anything else that I am missing to get the best out of what I can get?

Thank you very much in advance!
Alex

Answer: Hello again alslonik !

Here I'll add my little pinch of salt and recommend you having a look on that great manual : https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5100563/ . Of course it was not tested on plant genome, but it help you to orient your choice concerning assembly strategy depending on the technology you used and on your sequencing depth.

I already saw that you wanted to give it a try to quickmerge so you may have already saw that manual. As I tried quickmerge myself with 2 different PacBio only assembly, I have to say I was really satisfied with the result of quickmerge concerning contiguity and completeness. As you did a Falcon assembly, you can try merging a Falcon assembly and a Canu assembly, it may give some improvement as well, if you have the time of trying that of course !

Cheers,

Roxane
Question: Hi Biostar,

Apologies for this basic question but I can not find a simple answer

I want to sequence several samples of mRNA paired end 2x150 and for each sample I want 100 million reads (100 million for forward and 100 million for reverse)

the supplier proposes to sequence 3 samples per lane of flow cell of the Hiseq 4000 technology.
So what I understood 3 samples x 200 million reads (100 million for forward and 100 million for reverse) = 600 million reads per lane

but when I looked for the performances Illumina Hiseq 4000 I found this technology can generate 312000000 Clusters per lane only. How can this technology generate 600 million reads per lane ??

Would it be possible to give me a simple answer Thank you in advance

Answer: Illumina traditionally counts two reads coming from a single cluster independently. `312000000` passing clusters give rise to `312000000 x 2 = Number of reads` per lane.

Please keep in mind this number if attainable with `excellent` quality genomic DNA libraries. Any deviation from this will result in lower yields.
Question: Hello Everyone,

I am trying to use the R package KaryotypeR to produce rainfall plots in R. I am having trouble reading my VCF data into R and getting it to work in KaryotypeR. I've got their tutorial to work (link here):
[https://bernatgel.github.io/karyoploter_tutorial//Examples/Rainfall/Rainfall.html][1]

I have a VCF with example rows like this (excluding the many ##header lines for readability:

##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
##SAMPLE=
##SAMPLE=
##FILTER== 25 base quality">
##INFO=
##INFO=
##vcfProcessLog_20180918.1=,InputVCFSource=,InputVCFVer=<1.7.3>,InputVCFParam=>
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	NORMAL	TUMOUR
contig01000014.1	1128	165247d6-b4df-11e8-bec4-89a28c540d42	T	C	.	PASS	DP=89;MP=8.8e-01;GP=5.0e-11;TG=TT/CTTTT;TP=8.8e-01;SG=TT/TTTTT;SP=1.2e-01;CLPM=0.00;ASMD=136.00	GT:FAZ:FCZ:FGZ:FTZ:RAZ:RCZ:RGZ:RTZ:PM	0|0:0:0:0:11:0:0:1:27:0.0e+00	0|1:0:0:0:19:0:5:0:26:1.0e-01
contig01000543.1	276	16524d8a-b4df-11e8-bec4-89a28c540d42	T	C	.	PASS	DP=67;MP=8.7e-01;GP=1.3e-01;TG=TT/CTTTT;TP=8.7e-01;SG=CT/CTTTT;SP=1.2e-01;CLPM=0.00;ASMD=57.50	GT:FAZ:FCZ:FGZ:FTZ:RAZ:RCZ:RGZ:RTZ:PM	0|0:0:0:0:3:0:1:0:12:6.2e-02	0|1:0:3:1:21:0:3:0:23:1.2e-01
contig01000692.1	2433	16524ee8-b4df-11e8-bec4-89a28c540d42	A	T	.	PASS	DP=74;MP=9.8e-01;GP=2.0e-02;TG=AA/AAAAT;TP=8.8e-01;SG=AA/AAATT;SP=9.6e-02;CLPM=9.00;ASMD=101.00	GT:FAZ:FCZ:FGZ:FTZ:RAZ:RCZ:RGZ:RTZ:PM	0|0:19:0:0:3:25:0:0:0:6.4e-02	0|1:11:0:0:5:11:0:0:0:1.9e-01
contig01000694.1	830	16525032-b4df-11e8-bec4-89a28c540d42	C	T	.	PASS	DP=114;MP=1.0e+00;GP=5.6e-14;TG=CC/CCCCT;TP=8.6e-01;SG=CC/CCCTT;SP=1.4e-01;CLPM=0.00;ASMD=111.00	GT:FAZ:FCZ:FGZ:FTZ:RAZ:RCZ:RGZ:RTZ:PM	0|0:0:26:0:0:0:42:0:2:2.9e-02	0|1:0:8:0:1:0:26:0:9:2.3e-01
contig01000694.1	872	16525172-b4df-11e8-bec4-89a28c540d42	G	T	.	PASS	DP=89;MP=1.0e+00;GP=2.9e-07;TG=GG/GGGTT;TP=7.4e-01;SG=GG/GGGGT;SP=2.4e-01;CLPM=0.00;ASMD=115.50	GT:FAZ:FCZ:FGZ:FTZ:RAZ:RCZ:RGZ:RTZ:PM	0|0:0:0:23:0:0:0:34:3:5.0e-02	0|1:0:0:9:0:0:0:11:9:3.1e-01

I would like to convert this into a GRanges format such that I can subset the 'TUMOUR' data (excluding the 'NORMAL') and plot the variants as a rainfall plot with the KaryotypeR package using these commands from the tutorial (where sm.gr is a GRanges object of the data):

library(karyoploteR)
kp <- plotKaryotype(plot.type=4)
kpPlotRainfall(kp, data = sm.gr)

Could someone please point me towards a good way to read in and convert a VCF to a GRanges format suitable for this task?

Thanks for your help.

[1]: https://bernatgel.github.io/karyoploter_tutorial//Examples/Rainfall/Rainfall.html

Answer: library(karyoploteR)
tmp.vcf<-readLines("your_file.vcf")
tmp.vcf.data<-read.table("your_file.vcf")
tmp.vcf<-tmp.vcf[-(grep("#CHROM",tmp.vcf)+1):-(length(tmp.vcf))]
vcf.names<-unlist(strsplit(tmp.vcf[length(tmp.vcf)],"\t"))
names(tmp.vcf.data)<-vcf.names
colnames(tmp.vcf.data)[1] <- "Chr"
colnames(tmp.vcf.data)[2] <- "Start"
tmp.vcf.data = cbind(tmp.vcf.data,End=rep(tmp.vcf.data$Start+1))
tmp.vcf.data <- tmp.vcf.data[, c(1,2,12, 3:11)]
toGRanges(tmp.vcf.data)

Edit : Good luck to generate a custum genome based on your contigs in KaryoploteR !
Question: A lot of papers I've read have used [the qvalue package from Bioconductor][1] to do multiple testing correction. There is a function in the basic R stats package called [p.adjust][2] which seems to do the same thing. What is the advantage of the qvalue package, and why is it so frequently used instead of the built in R version?

[1]: http://www.bioconductor.org/packages/release/bioc/html/qvalue.html
[2]: https://stat.ethz.ch/R-manual/R-patched/library/stats/html/p.adjust.html

Answer: `p.adjust()` and the qvalue package aren't actually doing exactly the same thing. They are doing quite similar things with similar ends in mind, but the algorithms are different, so they produce different results (Gordon Smyth, who wrote `p.adjust()`, has a short summary of the history [here][1] and Tim Triche gives a nice explanation of the differences in [the same thread][2]). Interestingly, the packages I use almost always use `p.adjust()` rather than the qvalue package, though if you read a paper that uses them the authors might actually say "q-value" rather than "BH adjusted p-value".

BTW, q-values are a bit more difficult to calculate and you actually need quite a few data points. In effect, you have to use the data at hand to estimate the expected null rate and then multiply that by the BH adjusted p-value. Consequently, `p.adjust()` ends up being more generally applicable, but in cases where q-values are appropriate they tend to give more power.

[1]: https://stat.ethz.ch/pipermail/bioconductor/2012-December/049902.html
[2]: https://stat.ethz.ch/pipermail/bioconductor/attachments/20121219/00dc27b1/attachment.pl
Question: Hello,

If I have a set of orthologous proteins, it is assumed that these will have a shared function due to the maintenance of evolutionary pressure following speciation (as I understand).

If I have a set of proteins that fall within a single PFAM, is the assumption of shared functionality as strong as it would be for the ortholog group? I am guessing not, as the second paragraph of the reference's introduction explicitly says they do not always share function. But then, no assumption is complete, so I imagine there will also be orthologs that do not share function.

I assume there is a degree of uncertainty regarding how high quality the group of orthologs is (e.g. if it is found by RSD or RBH or tree based methods). But with that in mind, would you say the assumption of shared functionality is stronger for the PFAM group or for the orthologous group?

Thanks!

Reference:
Nucleic Acids Research, Volume 47, Issue D1, 08 January 2019, Pages D427–D432, https://doi.org/10.1093/nar/gky995

Edit - for anyone interested there's some interesting discussion here - https://www.biostars.org/p/128/#128

Answer: I would say stronger for the orthologous groups than in the PFAM groups. The pfam works only on conserved domains, while orthogroups (usually) take the whole sequence into account.

Moreover, it is not because you share a certain domain that you then also share a common function. Most eukaryotic proteins have a multi domain structure and thus the true biological function is a combination of those domains rather than 1 specific domain (though that also exists in eukaryotes). Some domain are very common (kinase domain, DNA binding, ... ) and others are quite specific.

If solely the absence or presence of a domain would determine function (and/or orthology) the problem of orthology defining would have been solved already (which it is not unfortunately)
Question: I am trying to plot a Heatmap. My data is as follows.

GO	s_3d	s_6d	s_17d
peptide metabolic process	4.71E-29	2.84828E-20	2.70719E-24
organonitrogen compound biosynthetic process	3.38158E-22	3.26674E-15	6.3903E-22
positive regulation of gene expression	0.057421474
negative regulation of RNA metabolic process	0.05933718
G-protein coupled receptor signaling pathway	0.061777665
positive regulation of macromolecule biosynthetic process	0.061906369
negative regulation of nitrogen compound metabolic process	0.065645467

Basically the data has 4 columns and a LOT of NA values.
I am using the function Heatmap.2 which usually has an option called na.color which can help me color the NA values to a distinct color

Since I have too many NA values, the heatmap function gives me an error

Error in hclustfun(distr) : NA/NaN/Inf in foreign function call (arg 11)

Therefore I tried to replace all NAs with zero and then give that a separate color. However the heatmap gives me a gradient even though I give the color separate. What I want is all 0's to be White and Everything which is not 0 to have a gradient from yellow to red. I do not want a gradient between white to yellow. How do I do this?

This is my code

rm(list = ls())
cat("\014")

if (!require("RColorBrewer")) {
install.packages("RColorBrewer", dependencies = TRUE)
library(RColorBrewer)
}

file <- "C:/Users/niran/Downloads/go_cluster.txt" #input
data <- read.table(file,header=T, sep = "\t",stringsAsFactors = F)
row.names(data) <- data[,1]
data <- data[-1]

# Since data is spread too far apart we can do a log transform to reduce the gaps

data[is.na(data)] <- 0 # Find possible solution for this
data <- -log10(data + 0.000000000001)

# creates a own color palette from yellow to red
color.palette <- colorRampPalette(c("#F8F8F8","yellow", "orange", "red"))(n=600)

col_breaks = c(seq(0,0.1,length=1), # for white
seq(0.1,2,length=100),
seq(2.01,4,length=100),
seq(4.01,6,length=100),
seq(6.01,8,length=100),
seq(8.01,10,length=100),
seq(10.01,12,length=100))

Heatmap <- heatmap.2(as.matrix(data),
main = "Heatmap of p-values for GO",
density.info = "none",
trace = "none",
margins = c(5,28),
key.xlab = "log10 Values",
cexRow = 1,
cexCol = 1.5,
keysize=0.75,
col = color.palette, # use on color palette defined earlier
breaks=col_breaks, # enable color transition at specified limits
dendrogram ="both")

dev.off()

Answer: Hi Seigfried,

you give the separate colors, but **ColorRamPalette** still interpolates them to create a smooth gradient from white to yellow to orange to red. If I understood you correctly, you would like first create a smooth gradient from yellow to orange to red and **then** prepend "white". Thus your code should contain something like:

color.palette <- c("#F8F8F8", colorRampPalette(c("yellow", "orange", "red"))(n=599))
Question: Hi,

We have run a pilot RNA-Seq study with one sample per condition, this is just a test run. I understand there is no valid statistical test in this case, however just curious to obtain differential expression through edgeR package in R assuming dispersion = 0.4 for the human data. I have a normal (baseline) sample followed by 5 different samples. When I run the edgeR package, I want to indicate my normal sample (as baseline), however, I am unsure what sample here is taken as a baseline for calculation of fold change.

**FC calculation**

> FC = Normal/Test_1 (OR)
> FC = Test_1/Normal

**Samples**

> Normal (baseline) = Normal
> Test (Treated) = Test_1
>

**Data**

dput(df_data)
structure(list(Normal = c(0L, 184L, 60L, 0L, 7L, 0L, 87L, 0L,
0L, 21L, 193L, 29L, 0L, 0L, 3L, 50L), Test_1 = c(0, 140.5, 64,
0, 4, 0, 83, 0, 1, 51.5, 199, 25, 0, 0, 5, 62)), class = "data.frame", row.names = c("Gene1",
"Gene2", "Gene3", "Gene4", "Gene5", "Gene6", "Gene7", "Gene8",
"Gene9", "Gene10", "Gene11", "Gene12", "Gene13", "Gene14", "Gene15",
"Gene16"))

dput(df_metadata)
structure(list(SampleID = c("xxxx1", "xxxx2"), CoreLabID = c("Normal",
"Test_1")), class = "data.frame", row.names = c("Normal", "Test_1"
))

**Here is the code that I am running**

bcv <- 0.4
y <- DGEList(counts=df_data, group=df_metadata$CoreLabID)
et <- exactTest(y, dispersion=bcv^2)

> View(et$table)
> structure(list(logFC = c(0, -0.67280976110796, -0.190706878123648,
> 0, -1.06592239047733, 0), logCPM = c(0.456451013758882, 6.84518828528986,
> 5.46338499556895, 0.456451013758882, 2.37389406911164, 0.456451013758882
> ), PValue = c(1, 0.433579402199822, 0.851984371429117, 1, 0.542580328250669,
> 1)), row.names = c("Gene1", "Gene2", "Gene3", "Gene4", "Gene5",
> "Gene6"), class = "data.frame")
>
> View(et$comparison)
> c("Test_1", "Normal")

Thank you,

Toufiq

Answer: edgeR has compared `Test_1` to `Normal` so, yes, the normal sample has been taken as the baseline. `et$comparison` records the numerator and denominator for the fold-change respectively.

You can find this out for yourself by reading the help page for `exactTest`, `help(exactTest)`. See the documentation for the `pair` argument.

If you use the standard edgeR pipeline
```
topTags(et)
```
then the results will be presented in a easy to read form and the output will tell you which sample has been compared to which.