cannin/biostars_qa huggingface dataset:
Question: Hi, I would perform a GO enrichment analysis for a set of genes, for exemple if I have 3 or 4 genes how to get the GO terms using topGO? in this [link][1] how to put my genes in 'allgenes' of the object topGOdata? Thank you

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

Answer: **NB - another answer here: https://www.biostars.org/p/325421/**

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

The input to topGO is a named list of genes and P-values, like this:

head(genes)
ANXA2 DHCR24 GBE1 GLA PRDX1 PSMD14
0 0 0 0 0 0
tail(genes)
CTTNBP2NL METTL7B AKR1C3 SLC5A11 TIMP4 PPP2R2C
7.395843e-07 1.496372e-06 4.328879e-05 1.612348e-04 2.620590e-04 9.728924e-04

range(genes)
[1] 0.0000000000 0.0009728924

The p-values are used to rank the genes, which is important when using the Kolmogorov-Smirnov test.

You can then run topGO like this (here with GO BP as database against which enrichment is performed).

require(topGO)
require(org.Hs.eg.db)

selection <- function(allScore){ return(allScore < 0.05)} # function that returns TRUE/FALSE for p-values<0.05
allGO2genes <- annFUN.org(whichOnto="BP", feasibleGenes=NULL, mapping="org.Hs.eg.db", ID="symbol")
GOdata <- new("topGOdata",
ontology="BP",
allGenes=genes,
annot=annFUN.GO2genes,
GO2genes=allGO2genes,
geneSel=selection,
nodeSize=10)

In order to make use of the rank information, use Kolmogorov-Smirnov (K-S) test:

results.ks <- runTest(GOdata, algorithm="classic", statistic="ks")
goEnrichment <- GenTable(GOdata, KS=results.ks, orderBy="KS", topNodes=20)
goEnrichment$KS <- as.numeric(goEnrichment$KS)
goEnrichment <- goEnrichment[goEnrichment$KS<0.05,]
goEnrichment <- goEnrichment[,c("GO.ID","Term","KS")]
goEnrichment$Term <- gsub(" [a-z]*\\.\\.\\.$", "", goEnrichment$Term)
goEnrichment$Term <- gsub("\\.\\.\\.$", "", goEnrichment$Term)
goEnrichment$Term <- paste(goEnrichment$GO.ID, goEnrichment$Term, sep=", ")
goEnrichment$Term <- factor(goEnrichment$Term, levels=rev(goEnrichment$Term))

Plot the results (the enrichment score is just negative log (base 10) of the enrichment P-value):

require(ggplot2)
ggplot(goEnrichment, aes(x=Term, y=-log10(KS))) +
stat_summary(geom = "bar", fun.y = mean, position = "dodge") +
xlab("Biological process") +
ylab("Enrichment") +
ggtitle("Title") +
scale_y_continuous(breaks = round(seq(0, max(-log10(goEnrichment$KS)), by = 2), 1)) +
theme_bw(base_size=24) +
theme(
legend.position='none',
legend.background=element_rect(),
plot.title=element_text(angle=0, size=24, face="bold", vjust=1),
axis.text.x=element_text(angle=0, size=18, face="bold", hjust=1.10),
axis.text.y=element_text(angle=0, size=18, face="bold", vjust=0.5),
axis.title=element_text(size=24, face="bold"),
legend.key=element_blank(), #removes the border
legend.key.size=unit(1, "cm"), #Sets overall area/size of the legend
legend.text=element_text(size=18), #Text size
title=element_text(size=18)) +
guides(colour=guide_legend(override.aes=list(size=2.5))) +
coord_flip()
Question: Hi !

Anyone know some linux commands that can be used to parse through the output coords table of [MUMmer][1]? It has the following output format:

[S1] [E1] | [S2] [E2] | [LEN 1] [LEN 2] | [% IDY] | [TAGS]
=====================================================================================
28536 30543 | 2008 1 | 2008 2008 | 89.71 | 3583	1233
6306 6469 | 972 808 | 164 165 | 92.73 | 3585	1364
6653 7101 | 455 4 | 449 452 | 84.31 | 3585	1364

All I need is the average of the % IDY column but have not been able to get this via linux commands. I also tried converting to an easier format such as txt but that also did not work.

Any ideas?
Thank you!

[1]: http://mummer.sourceforge.net/manual/#nucmer

Answer: as it looks your output is a set of fixed-length columns, you'd better try with something like:

cut -c 63-70
Question: When preforming a high-throughput RNA sequencing (with human samples), how many of the 3 billion base pairs of DNA will get covered after alignment and quality control?

In other words - how much of the genome can theoretically be reproduced from the RNA-seq?

I'm looking for even just a rough estimate, but if it helps, the samples that I'm interested in are from human brains, expressing ~16,000 genes, of which ~13,000 are protein coding.

I couldn't find an answer by googling, and will appreciate any help.

Answer: I have looked at a high-quality brain sample that I had worked on.

samtools depth -a data.bam | awk ' $3>0 { count += 1 } END { print (count/NR) }'

It appears that the genome is covered at about `7.5%` rate (command took 20 minutes to compute that)

It is true that the larger fraction of a genome may be observed transcribed into RNA, but those would not be picked up by a typical RNA-Seq experiment.
Question: Hi guys
I just can not figure out how to use the selectLab option in the biplot of PCAtools.
I tried creating a logical vector (withF for every sample I don't want to have a lab for) as well as a vector with all the sample_ids that should have a label. But nothing works nd I can't wrap my head around how to make it work.

Thanks a lot for any comments!
Sebastian

Answer: [ATpoint][1] and [Kevin Blighe][2] brought the solution!

selectLab expects a character vector with the exact labels we are handing to lab, not the corresponding row names.
If your metadata are factors, labs need as.character(), otherwise it prints the factor levels only.

biplot(data_PCA,
colby = "date_processed",
legendPosition = "bottom",
lab = as.character(data_PCA$metadata$genotype),
selectLab = c("Genotype1", "Genotype2"), #select here from $metadata$genotype!
pointSize = 6,
title = "PCA title",
caption = 'There is a clear seperation by date_processed')

Thanks a bunch guys, without you I couldn't have solved this!

[1]: https://www.biostars.org/u/25721/
[2]: https://www.biostars.org/u/41557/
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: I second ATpoint -- sounds like you need to familiarize yourself with how a matrix is interpreted in R, which is what `exprs(ALL)` returns.

A matrix is a vector with assigned dimensions, here these dimensions happen to be 2 (rows, columns).
So, when you call `median(matrix)`, it will simply take the median of all values within the matrix, ignoring the fact that it has rows and columns. If you want patient- and gene-specific values, you will have to specify that you want the values to be calculated on the rows (for genes) or columns (for patients), as ATpoint has shown.

Here's a toy example to illustrate the point:

```
> test_mat <- matrix(c(rep(1,5), 1:5,1:5), nrow = 5)
> test_mat
[,1] [,2] [,3]
[1,] 1 1 1
[2,] 1 2 2
[3,] 1 3 3
[4,] 1 4 4
[5,] 1 5 5

# median of all values
> median(test_mat)
[1] 2

# median for the values of the second column
> median(test_mat[,2])
[1] 3

# median for the values of the fifth row
> median(test_mat[5,])
[1] 5

# median calculated for every row
> apply(test_mat, 1, median)
[1] 1 2 3 4 5

# median calculated for every column
> apply(test_mat, 2, median)
[1] 1 3 3
```
Question: I'm about to look for piRNAs in some miRNAseq data. What are the most well established programs or pipelines for detecting piRNAs in human small RNA seq data?

Answer: I would not say there is one that stands out, but there are a couple which I would recommend:

- https://github.com/bowhan/piPipes
- https://www.smallrnagroup.uni-mainz.de/software.html
Question: how can i remove duplicated variants from vcf file? i googled and searched in biostars history but i did not fond any way to do it.

Answer: > . in fact i intend to remove variants that are same in scoffold id and pos and keep one of them.

I strongly suggest you also use the REF information...

sort on CHROM/POS/REF.
using awk create a KEY=CHROM\tPOS\REF, print the line if the key wasn't found previously

LC_ALL=C sort -t $'\t' -k1,1 -k2,2n -k4,4 input.vcf |\
awk -F '\t' '/^#/ {print;prev="";next;} {key=sprintf("%s\t%s\t%s",$1,$2,$4);if(key==prev) next;print;prev=key;}'

edit: added 'next; ' for VCF header.
Question: I am going to start GWAS analysis for human genome. Any recommendation which tools should I use? Also, it would be great if you could please send me a link to tutorial with mock data so I can practice first before working on the real data.

I am new to GWAS analysis so please help me. Many thanks in advance.

Answer: Read about [plink][1] and search for "plink tutorial GWAS". Also, related biostars post:

https://www.biostars.org/p/45440/
https://www.biostars.org/p/249149/

There are many R packages, too:

https://www.ncbi.nlm.nih.gov/pubmed/28132152
https://bioconductor.org/packages/release/bioc/html/GWASTools.html
https://github.com/JoniColeman/gwas_scripts
https://poissonisfish.wordpress.com/2017/10/09/genome-wide-association-studies-in-r/

[1]: https://www.cog-genomics.org/plink2
Question: Hi,

Is there a way to get a list of the start and end base pair positions of all known enhancers in the human genome (hg19) ?

Thanks in advance!

Answer: There may be other answers but one option would be table browser from UCSC (http://genome.ucsc.edu --> Tools menu --> Table Browser).
You would want to select `group: regulation table:vista enhancers` and then export the data in BED format.
Data comes from: http://enhancer.lbl.gov/
Question: Hello,
I am trying DELLY v0.8.2 for the first time, using ~x40 coverage data on a plant genome. The bam file was obtained by running bwa mem, then sorting and indexing the result.
I am using default params for DELLY. My run results in a bcf file with no data (i.e. header is present, but no variants detected) and prints the following message to stdout:

Sample:SRR1945819_map_to_pan.sort,ReadSize=51,Median=407,MAD=11,UniqueDiscordantPairs=230596

I find it very unlikely that there are no SVs between my sample and the reference, so I guess something is wrong with my data or the way I use DELLY.
Any ideas on what might be going wrong?
Thanks!

Answer: OK, this is a bit embarrassing, but turns out the problem was not with DELLY (or Manta), but rather with bcftools. Turns out I was using a very old version (0.1.19) that apparently can't read the new bcf. Once upgraded to 1.10.2, I am now able to read the output bcf, which is not empty at all.
Question: Hi ,

I am Aswathi currently working in V1-V3 amplicon sequencing and analysis. I am doing my sequencing in illumina platform, and as it is not possible to sequence ~515 bp length V1-V3, I did a trial to make two libraries of the V1-V3 amplicon.

One PE library of 2X250 bp covering outer region of the amplicon and second one covering inner region of the same amplicon. Therefore, inner library R1 and R2 reads are overlapping and outer R1 and R2 lies far away.
Then I tried to make a overlapped consensus amplicon sequences from these inner R1R2+Outer R1R2 merging together. I have used clustalO for the same. But I am able to produce very less (~5000) amplicon sequences using this methodology.

Can anyone suggest a better tool or approach for making amplicon from these two libraries?

Thank you.
Aswathi

Answer: If you search "paired end merge" you may get some information about software for this task.

FLASh worked well in my case. ( http://ccb.jhu.edu/software/FLASH/ ).

You could also try PEAR:
http://sco.h-its.org/exelixis/web/software/pear/

Well, these two tools only merge paired reads in one library. That means the IDs of each pair is same except /1 or /2 to distinguish the end. These tools may work for inner amplicon, but not for the outer amplicon.
These two tools are NOT for merging to libraries.
To merge two libraries, I would try:
1)	Merge inner amplicon with FLASh or PEAR to get Inner_merged set (Set1)
2)	Search Set1 with R1 and R2 of outer amplicon, respectively (using blast, blat or other tools)
3)	Write one script to get the best hit. R1R2 of outer amplicon, should hit the same fragment of Set1, and create one “faked name” for each pair.
4)	Merge R1_outer and Set1 to get Set2
5)	Merge Set2 with R2_outter.
Question: HI,

I need some help to start working with SAMtools. As a beginner with biostatistics background, it is my first time that I'm using Linux Ubuntu (I have installed it on my windows by making a machine on VMware for illustration purpose). Now I'm going to know how to work with SAMtools to carry out basic commands such as:

1- importing bam/sam files into samtools
2- Viewing the data
3- Sorting bam files
4. Calling variants

I have already searched on the internet and have found some tutorials such as http://www.htslib.org/doc/samtools.html , but I don't have any idea what those commands are!

could anybody please tell me where I can find a step by step tutorial for a beginner who even doesn't know how to read such commands?

If I am right, I should run samtools on terminal... by using the commands : "sudo apt-get install samtools" and then writing "samtools", I think I have run samtools but how should I write the commands in samtools on terminal?

I would also be grateful if someone tells me where I can find some low sized bam files to test samtools.

Thanks a lot!

Answer: It sounds like you are missing a lot of basic Unix knowledge that will make it difficult to progress without. If you are struggling with basic concepts such as terminal commands and their meaning, this will require you spend time and learn more about how to use a terminal in a Linux environment, and then spend some time reading the manual for the different programs you intend to use.

I'm going to go ahead and encourage you to think about using [Galaxy][1] for your project. It's a open, web-based platform for a variety of data analysis tasks and best of all features a GUI component for those that aren't very good or are rather new a terminal. [Here][2] seems to be a pretty good guide for variant calling using Galaxy, though admittedly i've never preformed this type of analysis so I cannot be entirely sure.

I also suggest that you somehow get a local machine with a copy of a Linux distro such as Ubuntu installed (perhaps ask your PI if this is possible), and not just a windows machine running a virtual machine. This will force you to get used to the Linux operating system and will make installing bioinformatic software such as samtools and eventually learning the terminal and the command line easier.

Here are a couple of links to get you started navigating the terminal and learning about the shell:

[Codeacademy][3]

[Command-Line Crash Course][4]

[Learning The Shell][5]

[Learn To Shell The Hard Way][6]

[1]: https://galaxyproject.org/
[2]: http://www.biomedcentral.com/content/supplementary/1756-0500-7-314-s1.pdf
[3]: https://www.codecademy.com/learn/learn-the-command-line
[4]: http://cli.learncodethehardway.org/book/
[5]: http://linuxcommand.org/learning_the_shell.php
[6]: http://learnpythonthehardway.org/book/appendix-a-cli/introduction.html
Question: Hello,
Can somebody tell me how to uncompress 1000 Genome vcf.gz files?
I am performing an RNA-editing analysis and would like to substract annotated SNPs/INDELs.
I have already done so using dbSNP data with bedtools instersect, but am still stuck with the 1000 Genome Project *.vcf.gz files.
I downloaded these for each chromosome and then concatenated them. These files are in a format that gunzip/gzip -d wont recognize.
I tried using this file unzipped in bedtools intersect but it wasn't reconized.
Many thanks,

Answer: try to use

bgzip -d

from the tabix package. But it's strange, gunzip should work.
Question: hi o all i have data like this

Ind	M1	M2	M3	M4	M5
P1	A/A	Unused	G/A	T/T	T/T
P2	T/T	A/A	A/A	A/A	G/G
1	T/A	A/A	G/A	T/T	G/G
2	T/T	A/A	G/A	T/T	T/G
3	T/T	A/A	G/A	T/T	T/G
4	T/T	A/A	G/A	A/T	G/G
5	T/A	A/A	G/A	A/T	T/G

First:i want to check whether P1 and P2 match each other for M1 or not, if they not matched(like P1=A/A,P2=T/T) that column should be selected and write in outfile, like this i want to do for remaining columns till M5 2nd step: Now in outputfile I want to match M1 value of 1st individual to P2, if it match then 1 otherwise it will be 0, M1 value is between them like (A/T ) then it will be H and so on for M4 and M5 3rd step: Then I will sum all the 1's and H'
finally my output will look like this

Ind	M1	M4	M5
P1
P2	SUM of 1's	SUM OF H
1	H	0	1	1	1
2	1	0	H	1	1
3	0	0	H	0	1
4	1	H	1	2	1
5	H	H	H	0	3

I am just beginner in R but till now i am doing this work in excel with sumif and countifs functions but now i can not do this work in excel because of column limitations. can any one help me to sort out this? any help in this regard will be highly appreciated. Here is my file for better understanding dropbox.com/s/ft943u7rf6v9j4b/row%20match.xlsx?dl=0
thanks in advance

Answer: Dear Friederike

Thank you very much for your kind and prompt reply to my request. I just modified your code little bit, till last step it is working fine but finally I am unable to format off_par. Now I am sharing my data file (https://www.dropbox.com/s/tysf6p3xdin9lop/biostars.xlsx?dl=0) and code I tried here (https://www.dropbox.com/s/l4ez5ht7i3o6xmd/Friederike%20Code.txt?dl=0) here for better understanding.
With kind regards,
Question: Hello, everyone. I searched in ENCODE for ChIP-Seq experiments and visualized them in UCSC browser. For instance, I searched for [TATA-box binding protein (TBP) assays][1] and went to [FOS' promoter][2] to examine its TBP peaks.

I kind-of get the idea of optimal Irreproducible Discovery Rate (IDR) for the peaks, but **I don't really understand the measure of "fold change over controls"** that UCSC gives you along the optimal IDR. Correct me if I'm wrong, but in the ChIP-Seq experiments' metadata (from ENCODE) it was stated that there was *no perturbation* (no treatment, no illness, nothing) in the samples, which makes me think those are the **controls**. If that's the case, then what exactly do they mean with "FC over controls", if the peaks are comming from "control" samples?

[1]: https://www.encodeproject.org/matrix/?type=Experiment&status=released&assay_slims=DNA%20binding&searchTerm=TBP
[2]: http://genome.ucsc.edu/cgi-bin/hgTracks?db=hg38&lastVirtModeType=default&lastVirtModeExtraState=&virtModeType=default&virtMode=0&nonVirtPosition=&position=chr14%3A75278096-75280364

Answer: Control most likely represents the chromatin input control that the peaks were called against, so the background levels of the assay.
If this is new to you I suggest you read about ChIP-seq and what kind of controls are commonly included.
Question: Hi, I have multiple sample bams which are split over two lanes.
e.g. for 2 samples:

flowcells: ["CCGB1ANXX-", "CCHNTANXX-"]
samples: ["3666-02-45-1_S2", "3666-03-45-1_S3"]
lanes: ["L005", "L006"]

I am trying to use snakemake to merge these bams. The problem I have is that it is not treating the sample ids separately in the merge rule but rather reading all sample bams and merging into one (for how ever many input samples there are so for two samples it is repeating the 'wrong' rule twice, for three samples, three times etc). How do I correct the rule so it only merges bams with the same sample id? Also, I have actually sorted these before this rule - should I sort after merging?
My code is:

rule mergeBams:
input:
bams = expand(intermediateDirectory + "{flowcell}{sample}_{lane}.unsorted.bam",flowcell = flowcells, sample = samples, lane=lanes)
output:
pileup = intermediateDirectory + "{sample}.merged.bam"
log:
"logs/{sample}.germlinePileup.log"
threads:
8
resources:
mem = "8GB",
time = "10:00:00"
container:
singularityContainers + "samtools-1.9.simg"
shell:
"samtools merge -@ {threads} {output.pileup} {input.bams}"

The std output from snakemake when running (it runs the rule twice):

samtools merge -@ 8 3666-03-45-1_S3.merged.bam CCGB1ANXX-3666-02-45-1_S2_L005.unsorted.bam CCGB1ANXX-3666-02-45-1_S2_L006.unsorted.bam CCGB1ANXX-3666-03-45-1_S3_L005.unsorted.bam CCGB1ANXX-3666-03-45-1_S3_L006.unsorted.bam CCHNTANXX-3666-02-45-1_S2_L005.unsorted.bam CCHNTANXX-3666-02-45-1_S2_L006.unsorted.bam CCHNTANXX-3666-03-45-1_S3_L005.unsorted.bam CCHNTANXX-3666-03-45-1_S3_L006.unsorted.bam

.....

samtools merge -@ 8 3666-02-45-1_S2.merged.bam CCGB1ANXX-3666-02-45-1_S2_L005.unsorted.bam CCGB1ANXX-3666-02-45-1_S2_L006.unsorted.bam /CCGB1ANXX-3666-03-45-1_S3_L005.unsorted.bam CCGB1ANXX-3666-03-45-1_S3_L006.unsorted.bam CCHNTANXX-3666-02-45-1_S2_L005.unsorted.bam CCHNTANXX-3666-02-45-1_S2_L006.unsorted.bam CCHNTANXX-3666-03-45-1_S3_L005.unsorted.bam CCHNTANXX-3666-03-45-1_S3_L006.unsorted.bam

Answer: If i understand correctly, you want to merge sub-samples with the same sample id to one bam file right?

If that's the case, you should use double bracket "{{}}" to let expand function keep the "{sample}" string instead of replacing them. Try the following input rule

input:
bams = expand(intermediateDirectory + "{flowcell}{{sample}}_{lane}.unsorted.bam",flowcell = flowcells, lane=lanes)
Question: Hi,

I'm a medical student with genetics, and I have come to understand some knowledge in bioinformatics would help me getting a research position in the field. I've been browsing the internet for some course, but all the ones I've found so far did not suit my needs. Does anyone know of an online bioinformatics course that:

1. is affordable
2. is open to an undergraduate medical student
3. provides a certificate
4. teaches you the most important things for genetic, genomic and molecular biology research

Answer: https://www.coursera.org/