cannin/biostars_qa huggingface dataset:
Question: Gene Ontology web specifically states that "a biological process is not equivalent to a pathway". They also state in their mission statement that "The mission of the GO Consortium is to develop an up-to-date, comprehensive, computational model of biological systems, from the molecular level to larger **pathways**, cellular and organism-level systems." and mention the term pathway in other sections.

I assume that GO Biological Processes are not equivalent to pathways but how are Biological Processes and pathway related? are biological processes part of pathways?

For example [Glycolytic process (GO:0006096)][1] and [Glucogenesis (GO:0006094)][2] terms seems to be pretty similar to the [Glycolysis / Gluconeogenesis pathway (KO00010)][3]

[1]: https://www.ebi.ac.uk/QuickGO/GTerm?id=GO:0006096
[2]: https://www.ebi.ac.uk/QuickGO/term/GO:0006094
[3]: https://www.genome.jp/dbget-bin/www_bget?ko00010

Answer: The word pathway used in the sentence you quote from the GO consortium doesn't have the same meaning as the word "pathway" as it refers to the KEGG term you link to.

The term "pathway" is common-usage word that is used in any number of different ways by any number of different groups of people/databases, where as the term "GO Biological Process" has a very defined meaning (defined by the GO Consortium). There are many different databases that claim to provide information about pathways e.g. KEGG, WikiPathways, IPA, Pathway Commons, etc. Each one probably has some subtly different idea about what they mean by the word pathway (I'm pretty sure several havn't ever really considered what makes something a pathway). Its important to realise that "pathway" is a human defined concept, and they probably don't ever correspond to a concrete, biologically isolatable reality.

That said there some things that I think most people would think applied to a "pathway", that other gene sets done. In a gene set defined by a shared GO term, nothing is implied about the relationship between genes within the set, other than the property of sharing some connection to the function implied in the go term. I think its reasonable, on the other hand, that for something to be called a "pathway", there must be some information about the ordering and relationship between the genes in the pathway i.e. A activates B, which represses C etc. These relationships don't have to be linear, but I think the word "pathway" implies that they are directed, and that pathways have inputs and outputs.
Question: Hi,

I am following the link: http://amos.sourceforge.net/wiki/index.php/Minimus2

for the merging of two contig files from velvet. It is written on the link: "......for merging one or two sequence sets (S1,S2)." here, do sequence sets mean contig files in fasta format?

I am facing the following error during merging:

command is:

./minimus2 all_velvet -D REFCOUNT=213048

error:

```
The log file is: all_velvet.runAmos.log
Doing step 10: Building AMOS bank & Dumping reads
Doing step 11
Doing step 12
Doing step 13
Doing step 20: Getting overlaps
Doing step 21
Command: /usr/bin/show-coords -H -c -l -o -r -I 94 all_velvet.delta | /home/kamran/Desktop/amos-3.1.0/bin/nucmerAnnotate | egrep 'BEGIN|END|CONTAIN|IDENTITY' > all_velvet.coords exited with status: 1
```

Help required

Regards,
Waqas

Answer: I know this post is OLD. Just answering anyway, in case anybody else has problems with this.

I got the same error message when testing minimus2 recently, and it turns out this was NOT because no overlaps would be detected (would have hoped for a more informative message on that, anyway).
Turns out, the standard error-output of minimus2 is not very helpful. You have to check the ".runAmos.log" file for actual info on what went wrong. (in your case that would be "all_velvet.runAmos.log".

In my case I found that minimus was expecting the tools `delta-filter` and `show-coords` at `/usr/local/bin`. However this may not be the case if you installed amos locally on a shared server, where you have no admin rights (as I did via bioconda).

Therefore you may have to manually edit the minimus2 script at lines `46` & `47` to point to the correct executables.
Solved the problem for me!
Question: Hi all
i need help in one silly question
Can any one help me understand the meaning of p value?
if we say **p value = 0.05** then what does that means.
what is confidence interval?
If we say **95% confidence** interval then what does that means.

I read certain article where it says p value is error estimation . If we are dealing with **trancriptome data** then in this data what could be error. meaning people select **pvalue < 0.05 as significance** what this significance means???

Thanks ..

Answer: When running a test you are running a risk of type I error (accepting null hypothesis when it is not true, false positives, alpha) and type II error (failing to reject the null hypothesis when in fact it is false, false negatives, beta). P value is the smallest probability (alpha) at which your test is still significant. P-value of 0.05 basically means you are running a 5% risk that your test will be significant by chance and this value is commonly accepted as sufficient for a single test. So when running a test and returning the P-value of e.g. 0.003 (which is smaller than 0.05) you accept the value as significant since the probability of running a type I error (false positives, basically your P value) is smaller than prior selected (0.05).

Confidence intervals give you a range of values that would be significant at a preselected confidence level (1 - selected P value * 100 --> 1 - 0.05 * 100 = 95%).
Question: Hi,

I have a Python script which takes as input a list of SNPs, then calls tabix to retrieve .VCF file and does some stuff with it.
The problem is that i need to manually add (using UCSC genome browser web*) the coordinates of the SNP for tabix to work properly.

Does any one know any simple python API / function / module to easy add this information without having to perform the search manually?

Thanks!

Answer: You could use an external database to find the coordinates of your rs## (not tested, please check the +0/+1 for the coordinates)

for R in rs25 rs26 rs27
do
mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -Dhg19 -N -e "select concat(chrom,':',chromStart+1,'-',chromEnd) from snp135 where name='$R';" |\
xargs tabix my.vcf.gz
done
Question: I understand the meaning of N50 contig length, but lately I've been coming across assembly statistics that report the N50 contig number. Can someone explain what is meant by N50 contig number and how it relates to assembly quality?

Answer: N50-length is the length of that last contig that puts you over the top.
N50-contig number (i.e. N50) is just the rank of that contig.

91 77 70 69 62 56 45 29 16 4[?]N50 vs N50 length[?]
Technically N50, as opposed to N50 length, refers to the ordinal of that last contig that pushes it over the brink - in this example 4 (since 69bp is the 4th largest contig). Unfortunately, a higher N50 implies the opposite of a longer N50 length. Some papers refer to N50 length as L50, while most have simply followed the lazy convention of dropping "length" off of "N50 length". I think it is important to include units with your N50 to minimize confusion.

http://jermdemo.blogspot.com/2008/11/calculating-n50-from-velvet-output.html
Question: In UCSC genome browser web server there are some chromosome which start with "chrUn", for example chrUn_KI270438v1. What are these chromosome exactly? It means Unknown chromosome? If yes, why unknown?

Answer: Quoting directly from the UCSC website:

*"ChrUn contains clone contigs that cannot be confidently placed on a specific chromosome. For the chrN_random and chrUn_random files, we essentially just concatenate together all the contigs into short pseudo-chromosomes. The coordinates of these are fairly arbitrary, although the relative positions of the coordinates are good within a contig."*

https://genome.ucsc.edu/FAQ/FAQdownloads.html
Question: Hi Everybody,

I'm trying to run a gene regulatory inference method.
In the first step I need to remove genes with no specific change in their expression values.
I don't know how to specify if a gene has "significant" change. What does it mean?
My second question:
Should I just remove these genes (not significantly changed) if I'm looking for their regulators?
Or a gene which has not changed significantly can't be also detected as a transcription factor?
In other words:
I thought maybe a gene is not changed significantly, so I am not be able to find those genes regulating this gene.
But what if I want to detect the influence of this gene on other genes? If it is not changed, is it not worthy to consider its effect on other genes? And gene regulatory network inference methods only rely on the changes to find regulatory interactions?

Thanks.

Answer: You may want to use a variance filter:
(in R)

data$vars = apply(data, 1, var) #Calculate variance across rows

data2 = data[data$vars >= quantile(data$vars, c(.80)), ] #Eliminate rows, here using 80% percentile

data2$vars <- NULL #Removing the variance column

Perhaps someone with more R experience can rewrite this to a more concise filtering :-)
Question: There are many R packages from where GO annotations can be retrieved, for example, GO.db, biomaRt, and org.Hs.eg.db. Are there any differences between the annotations obtained from these resources? Which one is the best in terms of update frequency and easy to use?

One case example adapted from: https://support.bioconductor.org/p/38420/

> library(GO.db)

> library(org.Hs.eg.db)

> res <- get("GO:0006913", revmap(org.Hs.egGO))

> res <- do.call('c', mget(res, org.Hs.egSYMBOL))

> sort(unique(res))

[1] "AAAS" "ANKRD54" "ANP32A" "CAMK1" "CDK5" "CITED1" "EIF5A" "FBXO22" "MYBBP1A" "NEUROD1" "NPM1" "NSRP1" "NUP205"
[14] "NUP98" "RGS14" "RSRC1" "SET" "UPF3A"
>
> library(biomaRt)

> ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")

> res <- getBM(c("hgnc_symbol"), filters = "go", values = "GO:0006913", mart = ensembl)

> sort(unique(res[,1]))

[1] "AAAS" "ANKRD54" "ANP32A" "ANP32D" "ANP32E" "CAMK1" "CDK5" "CITED1" "EIF5A" "FBXO22" "MYBBP1A" "NEUROD1" "NPM1"
[14] "NSRP1" "NUP155" "NUP205" "NUP35" "NUP54" "NUP58" "NUP98" "RAN" "RGS14" "RSRC1" "SET" "UPF3A"

> sessionInfo()

R version 3.4.1 (2017-06-30)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Sierra 10.12.6

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib

locale:
[1] C

attached base packages:
[1] parallel stats4 stats graphics grDevices utils datasets methods base

other attached packages:
[1] data.table_1.10.4-3 gap_1.1-20 biomaRt_2.32.1 GO.db_3.4.1 org.Hs.eg.db_3.4.1 AnnotationDbi_1.38.2
[7] IRanges_2.10.5 S4Vectors_0.14.7 Biobase_2.36.2 BiocGenerics_0.22.1 limma_3.32.10

loaded via a namespace (and not attached):
[1] Rcpp_0.12.14 bit_1.1-12 rlang_0.1.4 blob_1.1.0 plyr_1.8.4 tools_3.4.1 DBI_0.7 bit64_0.9-7
[9] digest_0.6.12 tibble_1.3.4 bitops_1.0-6 RCurl_1.95-4.8 memoise_1.1.0 RSQLite_2.0 compiler_3.4.1 XML_3.98-1.9
[17] pkgconfig_2.0.1

Answer: GO.db and org.Hs.eg.db are copies of the GO annotations. GO.db is updated every 6 months with each release of Bioconductor. org.Hs.eg.db is also updated at the same time and using GO.db. Following the release schedule of Bioconductor as explained in [Annotation Packages][1]:
> Annotation packages contain lightly or non-curated data from a public source and are updated with each Bioconductor release (every 6 months).

biomaRt connects to the server where the informations is stored, so it will be the most up to date. Ensembl is updated every 2-3 months according to [Emily][2]:

> BioMart is not updated every day. BioMart is updated with the rest of the Ensembl database, every 2-3 months.

If you want a stable release you can use either GO.db or org.Hs.eg.db, if you want the most up to date data every time you do an analysis you can use biomaRt.

[1]: https://bioconductor.org/developers/package-submission/#submission
[2]: https://bioinformatics.stackexchange.com/questions/8799/discordance-of-gene-ontology-data-between-mouse-genome-informatics-and-ensembl-s/8801?noredirect=1#comment11858_8801
Question: Hi there,

I want to get controlled access to dbGaP - in particular to the data that is linked here (section: "Download from dbGaP"): https://github.com/joed3/GTExV6PRareVariation

I am based at the University of Potsdam, Germany and hence have none of the accounts requested for application (eRA Commons, DUNS number, account in SAM and Grant portal, ...) as I am not enrolled in any project funded by the US government or similar. To me, it seems as if I cannot get access at all to the data. I have read all the FAQ's and manuals I could find at dbGaP and eRA.

Is there anyone (from a non-US institution) who has walked through the application process of dbGaP and can help me with what accounts I really need and can create at all (i.e. to get an eRA Commons account, I should be fine with the rest)?

Best,

Cindy

Answer: Yes, we access lots of dbGap data from the UK.

The first step is to apply for an eRA account. You will need the help of whomever is designated as the "signing officier" at your university. If your university does not have one, then your university will first need to register. See here:
https://era.nih.gov/commons-account-information.cfm

Once your university has an assigned signing officer and account administrator, then you (assuming you are a PI) ask the SO to create you an eRA account associated with the University. You can then submit a request to dbGap (for which you will require help from the SO and in addition the universities head of information security/head of IT).
Question: bgzip files are backward compatible with gzip, but I have issues when using bgzip compressed vcf files with snpeff (java) or perl scripts that uses IO::Uncompress::Gunzip (that I believe it uses zlib under the hood). In both cases the data is decompressed but truncated after a few hundred lines aprox. I could be totally wrong but I was wondering if zlib (or whatever gzip compatible library they are using) is getting confused with the bgzip bloks and only processing one or a few of them leaving the output incomplete.

perl code that does not work:

#!/usr/bin/env perl
use strict;
use IO::Uncompress::Gunzip qw(gunzip $GunzipError) ;

my $infile = shift;
my $infh = IO::Uncompress::Gunzip->new( $infile ) or die "IO::Uncompress::Gunzip failed: $GunzipError\n";
my $line_count = 0;
while (my $line=<$infh>){

$line_count++
}
print "total lines read = $line_count\n";

This gives 419 lines

$ perl /home/pmg/tmp/test_zlib-bgzip.pl 460112_TTAGGC_L005_L006_C3HVJACXX.sorted.rmdup.varsit.vcf.gz
total lines read = 419

but using open with gzip pipe works:

#!/usr/bin/env perl
use strict;

my $infile = shift;
open(my $infh , 'gzip -dc '.$infile.' |'); # I can use bgzip intead gzip
my $line_count = 0;
while (my $line=<$infh>){

$line_count++
}
print "total lines read = $line_count\n";

Gives the expected number of lines

$ perl /home/pmg/tmp/test_gzip-bgzip.pl 460112_TTAGGC_L005_L006_C3HVJACXX.sorted.rmdup.varsit.vcf.gz
total lines read = 652829

I googled about and I was unable to find quickly any relevant entry, but this is something that I am sure other people would have already faced. Do someone have a clue about why is this happening? I am using ubuntu 12.04.4 with perl 5.16

Answer: See the SAM spec about the bug/feature in the Java library - it only reads the first gzip block. There is also a related thread in BioStar. I guess Perl is the same. On Unix/Linux, I usually use:

open(FH, $ARGV[0] =~ /\.gz$/? "gzip -dc $ARGV[0] |" : $ARGV[0]) || die;
while (<FH>) {...}
Question: Hello All,

I have ben trying to run Htseq-count for quite a while now. Its containing throwing error as follows,
Please Install PySam to use the BAM_Reader Class (http://code.google.com/p/pysam/) Error occurred when reading the beginning of SAM/BAM file.
No module named pysam

And so as it is mentioned I have installed pysam and when I run Htseq-count it still complains as No module named pysam
I have installed pysam using pip install command line and as well as downloading the file.
When I insatll pysam using pip its show as it is there
Requirement already satisfied (use --upgrade to upgrade): pysam in /home/user/Tools/anaconda3/lib/python3.4/site-packages/pysam-0.9.0-py3.4-linux-x86_64.egg

However, I think the htseq-count doesn't consider anaconda version of pysam. And so as I think its not taking into account the downloaded and installed pysam package also.

Any suggestions would be great on this issue

Answer: htseq-count will be using your default python installation and therefore not find the pysam module in your anaconda installation. He have had the same issue.

Two options:

-Install pysam for your default python installation

-change the python interpreter in the htseq-count script. Find the script using 'which htseq-count', edit the interpreter at the top of the script to the path to your anaconda python.
Question: I came across a slide share about the new GRCh38 assembly and its ALT_LOCI assemblies.
(http://www.slideshare.net/vaschn/agbt2014schneider)

Question:
1. Are there any "ALT Loci aware" aligners out there ?

If so, then how are accommodating these ALT LOCI while aligning reads to the reference genome ? How do they chose between primary assembly and Alternate Loci ?

Some references:
This Biostar thread (Applying patches to GRCh assembly) talks about patches/Alts during alignment as well as about the aligner - srprism (ftp://ftp.ncbi.nlm.nih.gov/pub/agarwala/srprism)

Nice overview of ALT LOCI in GRCh38 (Will GRC38/HG20 be a multiple sequence reference genome?)

thanks in Advance,

Aniket

Answer: There are two 'ALT loci aware' (or multi-allelic reference) aware aligners that I know about.

srprism (ftp://ftp.ncbi.nlm.nih.gov/pub/agarwala/srprism)

ibwa (https://github.com/genome/ibwa)

Neither are published (AFAIK) but srprism is in preparation. In general, I think they both need to be tested more.
It is also unclear how other tools along the tool chain (variant calling) will respond- and they will likely also need modification to take advantage of the full assembly.
Question: Hi all.

I was running coverageBed with my exome bam file but it takes too long time and it hasn't returned output value.

Is it right that coverageBed takes much time to finish?

my bed file was very tiny so I thought that something might go in a wrong way.

below is my command and my bed file.

coverageBed -abam /DATA1/dmcb/TAR22_WES/TAR22S_sequence_readgroup_sort_dup_remove_realigned_recal.bam -b fgfr.bed -d > TAR22.cov.del.b36

10 123239372 123239472

Answer: You might try this:

samtools view -u -L fgfr.bed myaln.bam \
| coverageBed -abam - -b fgfr.bed -d > out.bed

Which means: First get the reads overlapping your bed file, then do the counting.

With 20GB of bam file I'm not surprised coverageBed takes more then 10 min.

(By the way, the title of your post is misleading as there is no error reported in your question)
Question: I have tumor age (nominal variable, divided into age brackets), gene expression (in RPKM) and sex (binary) for many tumors in people at different ages (not paired..ie. not time series of same tumors but different tumors at each age). I want to find the genes that are most differentially expressed with age, controlling for sex. What would be the best way to do this?

More specifically:

1. Do I need to transform the RPKM values to rank (as per [this thread][1])? The difference is that they use regularized regression.
2. Is a linear regression appropriate? [This paper][2] seems to suggest quantile regression is better (however there are other features of their algorithm, and their age is not ordinal) to linear regression which is used in much of the aging gene expression studies with microarray.
3. Because I am repurposing data and do not know batch would it be appropriate to run sva (or something similar) on it?

Thanks!

[1]: https://www.biostars.org/p/105937/
[2]: http://www.biomedcentral.com/1471-2164/10/S3/S16

Answer: If possible, get the raw data, RPKMs have a LOT of issues.

1. Depending on how complex you design becomes, it can become difficult to analyze rank transformed data.
2. Give linear regression a try first, there are a lot more tools for it. I'd only deal with quantile regression if the results were unsatisfactory.
3. Absolutely. You might also try just including the batches as a factor, though odds are good that SVA will work a bit better.
Question: I am trying to combine 4 fastq.gz files in to one big file using following command.

cat *.fastq.gz > File_Big.fastq.gz

However Every time I got following error

Error: object '.fastq.gz' not found

I made sure that I am working in correct directory.

Anyone have any idea where I am going wrong?

Thanks

Answer: You are using this command inside **R**.

Don't use R directly, using a linux shell. Like bash. https://www.tldp.org/LDP/Bash-Beginners-Guide/html/Bash-Beginners-Guide.html
Question: Hi all,

I have this doubt which I couldn't sort out by myself. I started working with VCF files since one month only, so it might be that it's just that I didn't stumble on the right piece of literature yet. My VCF lines in the file have an INFO field which is like this:

DP=20;VDB=0.000678823;SGB=-1.30291;MQ0F=0.5;AC=4;AN=4;DP4=0,0,19,0;MQ=3

Some of these Key=Value pairs are totally fine for me, the documentation is very clear. However, I can't find the description for all of them. For example: VDB, SGB, I don't know what they are.

Where can I find such information?

Answer: you should find this information in the VCF header, in a INFO field

grep '##INFO=