cannin/biostars_qa huggingface dataset:
Question: Hi, I've collected annotations of a large number of genes from GO, and propagated annotations back to the GO ancestors. Now I am trying to use these annotations for a general-purpose genomic functional graph (similar to genomic KG, etc. that are recently developed). The problem is that those top-level terms (e.g., on levels 0 and 1) really don't make sense. However, I don't really know at which level I should be cutting off the GO terms... And there are also papers cutting off terms based on the number of genes associated with such terms.

Can someone recommend some papers or approaches that contain details on essentially trimming GO terms w.r.t GO level or #genes associated (whether before or after propagation of GOAs), etc.?

Maybe I should ask @go

Answer: Suggestion from GO staff:

>I would suggest you use one of the GO subsets available http://geneontology.org/docs/download-ontology/

>The 'Generic GO subset' is one that should be broadly useful for your case, giving enough granularity to map genes to meanful classes. Organism-specific subsets are also available.
Question: My goal is to split a sequence at a specific site into two separate sequences. Searching for the site should be a bit fuzzy due to sequencing-pipeline (basecalling on MinION) error.

Example:

Assume a sequence as below. X, Y, Q and Z are sequence nucleotides not necessary for understanding the problem but are useful for demonstration purposes.

XXXXXXXXXXXXXXYYYYYACTCATAQQQQQQQQQZZZZZZZZZZZZZ
|-----|

I would like to find site `ACTCATA` (with fuzzy matching) and split the sequence into

XXXXXXXXXXXXXXYYYYY

and

QQQQQQQQQZZZZZZZZZZZZZ

with optionally discarding the matched sequence.

Bonus points if this is done on fastq files where data on quality of reads is also split into new strings.

This could probably be accomplished the pedestrian way in biopython but was wondering if I missed a tool that does what I describe above.

Answer: [seqkit](https://github.com/shenwei356/seqkit)

$ cat read.fq
@seq
XXXXXXXXXXXXXXYYYYYACTCATAQQQQQQQQQZZZZZZZZZZZZZ
+
GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG

$ cat read.fq | seqkit locate -p ACTCATA
seqID patternName pattern strand start end matched
seq ACTCATA ACTCATA + 20 26 ACTCATA

$ cat read.fq | seqkit subseq -r 1:19
@seq
XXXXXXXXXXXXXXYYYYY
+
GGGGGGGGGGGGGGGGGGG

$ cat read.fq | seqkit subseq -r 27:-1
@seq
QQQQQQQQQZZZZZZZZZZZZZ
+
GGGGGGGGGGGGGGGGGGGGGG
Question: How can I download Exons coordinates & exon Numbers associated with each gene for only **one transcript per Gene**?I know there are many transcript per Gene. have tried biomart & Table browser but it gives me result for all the transcript.

Answer: Thank you All. I think i have found solution.
https://groups.google.com/a/soe.ucsc.edu/forum/#!msg/genome/_6asF5KciPc/VTm611kJHAAJ
Question: Hi, everyone!
Is there anyone here familiar with or using the RTCGAToolbox package?
I found it to be very useful when importing, compiling and exporting data from the Firehose repository. However, I am simply unable to export miRNA data.
Here is what I'm doing:

> miOV <- getFirehoseData(dataset="OV", runDate = "20160128", gistic2Date = NULL, clinical = TRUE, miRNAArray = T, forceDownload = T)
And it works fine:

> miOV
OV FirehoseData objectStandard run date: 20160128
Analysis running date: N/A
Available data types:
clinical: A data frame of phenotype data, dim: 591 x 21
miRNAArray: A list of FirehosemRNAArray object(s), length: 1
To export data, use the 'getData' function.

Now, when I try to export it, it gives me the following message:

> miRNAData <- getData(OV_miRNA, "miRNAArray")

> Error in .getListData(getElement(object, type), platform) :
argument "platform" is missing, with no default

I checked what the acceptable arguments were, and it describes as:

> An index for data types that may come from multiple platforms (such as mRNAArray), for GISTIC data, one of the options: 'AllByGene' or 'ThresholdedByGene'

In this case, I checked the platform for [this specific set of data][1], which was performed on Agilent Human miRNA microarray Rel12.0.

I tried adding "Agilent", "Agilent Human miRNA", and other values, but it didn't work.

Would anybody have any idea on how to solve this issue?
Thanks!

[1]: http://gdac.broadinstitute.org/runs/stddata__latest/samples_report/OV.html

Answer: Case solved!
Although it still gives an error message after performing the getData function, it still exports the data.
Question: Hi guys

Is it possible to buy / license a local copy of the entire pubmed / medline with full-text?

Best,
Kasper

Answer: No. You can get ~1 million open access articles from ftp://ftp.ncbi.nlm.nih.gov/pub/pmc. For the other 23 million, you need to get access to full text from the publisher. Check the Genocoding project at http://text.soe.ucsc.edu/ since they are trying to mine sequences from all Pubmed and have a table on what publishers have granted access.
Question: Hi everyone,

I have a list of genes (n = 5726) for which I wanted to get the chrom, start and end. I used the following query:

library(biomaRt)
ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
lit.genes.ann <- getBM(attributes = c('chromosome_name','start_position','end_position','hgnc_symbol'), filters = 'hgnc_symbol', values = gene.names, mart = ensembl)

I got results for all but 79 genes back. This list was probably generated some time back and I am using the current ensembl build to get the data. I looked up one such gene EMR1 on Genecards and looks like it has been changed to [ADGRE1][1].

Is there a way to do a batch search & import to get current HGNC symbols for all the obsolete genes that I have? I am not looking for all the aliases but just one approved HGNC symbol per obsolete gene name. Some of the gene names are:

ACN9
AQPEP
BAI1
PAPL
EMR1

[1]: http://www.genecards.org/cgi-bin/carddisp.pl?gene=ADGRE1

Answer: Alright, I posted the question too soon I guess. For future reference - [Multi-symbol checker][1] worked.

[1]: http://www.genenames.org/cgi-bin/symbol_checker
Question: Hi there, I'm new to this analysis area, and I already look into many threads, but still don't understand the meaning of "mean normalized counts" in MA plot. Is it the mean normalized read counts between comparison groups? If so, what's the biological meaning of it, and also with fold change ? Thanks!

Answer: It indicates the average expression between the groups you compared, which is simply the mean of the normalized counts of the samples used in the comparison group. The "meaning" is that genes with higher counts are more trustworthy (=have higher statistical power) than low count genes. Biologically this is because either the gene is lowly expressed and/or short and therefore intrinsically gets lower counts than a longer one at equal expression level.
Question: How to do SNP selection in large VCF/BCF/GenotypeTables in R like you can using BCFTools,CYVCF2 (fast python VCF/BCF parser) or Excel?

Support for at least the following type of filter criteria is required:

- Filtering on Region of Interest
- Filtering on Samples of interest
- Filtering on Variant Quality (QUAL)
- Filtering on Genotype call rate (INFO/AN / (number of samples*2))
- Number of alternative alleles (INFO/AC)
- Filtering on alternative allele frequency (INFO/AF)
- Samples required to HOM_REF,HET,HOM_ALT
- Samples to have contrasting genotypes

Preferred supported filtering criteria:

- Filtering on sample specific or minimal/maximum/average genotype depth
- Filtering on sample specific or minimal/maximum/average genotype quality
- Distance to nearest neighbor variants upstream and downstream

Depending on filter support for genotype depth and genotype quality these VCF/BCF/R specific format files can be quite big.

So compact (small data size) R specific format and multi-threaded filtering might also be useful/required when loading everything to memory on a big machine.

Or streaming parsing of the VCF/BCF file (like BCFtools and CYVCF2 do) .

What would be the best way and library to do SNP selection using the above filter criteria in R?

Is there maybe something like CYVCF2 (python wrapper around HTSLIB VCF parser) for R?

Or maybe a BCFTools like something in R?

Or some other good R library for doing this?

The targeted downstream users are R(studio) users who would like to have this functionality in the environment that they are used to and already use for related tasks.

The VCF can optionally be (pre) converted to TAB delimited format with bcftools query. Or to any specific R format.

A high level / coarse filter on region of interest and samples of interest filter might already have been done.
It's often (but not always) more about iterative / dynamic fine grained filtering final selection of SNPs.

Answer: Have a look at the VariantAnnotation package in Bioconductor.

https://bioconductor.org/packages/release/bioc/html/VariantAnnotation.html

And the vignette on filtering VCF files.

https://bioconductor.org/packages/release/bioc/vignettes/VariantAnnotation/inst/doc/filterVcf.pdf

Note that this may not be the fastest option due to the parsing into VCF object in R, but it can be a useful adjunct to "system"
calls (remember, you can wrap command-line utilities in R functions to make them "appear" to R users as R functions).
Question: Hello,

I have a genome with homologous pairs of chromosome. I am using hisat2 to map RNA-seq reads to genome. I want to do variant calling with GATK, but before I want to extract the reads that mapped uniquely from BAM file(not SAM). Could you suggest how to do that? From other threads seems that flag NH:i:1 means mapped once, but how can I grep it from bam file? and is it a right way?

HISEQ1:128:C0HUMACXX:6:2210:15441:67674	163	chr1A	1101543	60	59M	=	1104438	626	CAATCTTCTCGTGGTTTATTGGCGAATGGAGTTTCAGAAGATGCTCAAAATGGTTTCCA	@BCFFDDEFHDDFGEEEDHEHCBHFD<3D3C*?D:B?:B?B*98D?4?BB=CG)=FFHA	AS:i:0	ZS:i:0	XN:i:0	XM:i:0	XO:i:0	XG:i:0	NM:i:0	MD:Z:59	YS:i:0	YT:Z:CP	NH:i:1
HISEQ1:130:C0G55ACXX:1:1203:2254:30173	163	chr1A	1101543	60	101M	=	1101599	155	CAATCTTCTCGTGGTTTATTGGCGAATGGAGTTTCAGAAGATGCTCAAAATGGTTTCCATGATGCTGATCGAACTGTCCGCAGAGGTGAGGGGCCATCTAA	CCCFFFFFHHHHHJGIJJJJJJJJJJJJJGIJIJJIJIJJJJJJJJJJJJJJJIIJJJJJJJJJJJJHHHHHFFFFEEEDDDDDDD@BDDDDDDDDDDDDD	AS:i:0	ZS:i:-15	XN:i:0	XM:i:0	XO:i:0	XG:i:0	NM:i:0	MD:Z:101	YS:i:-2	YT:Z:CP	NH:i:1

Answer: `samtools view -b -q 5 -o output.bam alignments.bam`

Having said that, there's no reason to bother doing this. GATK is smart enough to understand that low MAPQ alignments should be ignored.
Question: Have anyone noticed the BioMart servers were down

I tried the UK, usesat, uswest and Asia. None of them gave response.

Sent a ticket to Helpdesk in Ensembl.

---

> listMarts(host='www.ensembl.org')
Error in listMarts(host = "www.ensembl.org") :
Unexpected format to the list of available marts.
Please check the following URL manually, and try ?listMarts for advice.
http://www.ensembl.org:80/biomart/martservice?type=registry&requestid=biomaRt

> listMarts(host='useast.ensembl.org')
Request to BioMart web service failed.
The BioMart web service you're accessing may be down.
Check the following URL and see if this website is available:
http://useast.ensembl.org:80/biomart/martservice?type=registry&requestid=biomaRt
Error in if (!grepl(x = registry, pattern = "^\n*")) { :
argument is of length zero

> listMarts(host='uswest.ensembl.org')
Request to BioMart web service failed.
The BioMart web service you're accessing may be down.
Check the following URL and see if this website is available:
http://uswest.ensembl.org:80/biomart/martservice?type=registry&requestid=biomaRt
Error in if (!grepl(x = registry, pattern = "^\n*")) { :
argument is of length zero

Answer: Hi,

I'm sorry that you've had a problem while using BioMart. BioMart experienced a heavy load yesterday's evening due to a number of large jobs monopolising our resources. I can confirm it's up and running now. Apologies for the inconvenience and thank you for your patience.

Best wishes,
Michal
Ensembl Helpdesk
Question: Hi all,

I am using clusterProfiler to perform KEGG enrichment. I have a list of gene Symbol (Rat genes). First I need to translate the symbols to EntrezID.

I installed the Genome wide annotation for Rat:

source("https://bioconductor.org/biocLite.R")
biocLite("org.Rn.eg.db")

Then I tried to run the command:

x <- c("GPX3", "GLRX", "LBP", "CRYAB", "DEFB1", "HCLS1", "SOD2", "HSPA2",
"ORM1", "IGFBP1", "PTHLH", "GPC3", "IGFBP3","TOB1", "MITF", "NDRG1",
"NR1H4", "FGFR3", "PVR", "IL6", "PTPRM", "ERBB2", "NID2", "LAMB1",
"COMP", "PLS3", "MCAM", "SPP1", "LAMC1", "COL4A2", "COL4A1", "MYOC",
"ANXA4", "TFPI2", "CST6", "SLPI", "TIMP2", "CPM", "GGT1", "NNMT",
"MAL", "EEF1A2", "HGD", "TCN2", "CDA", "PCCA", "CRYM", "PDXK",
"STC1", "WARS", "HMOX1", "FXYD2", "RBP4", "SLC6A12", "KDELR3", "ITM2B")

eg = bitr(x, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Rn.eg.db")

Then I got the warning message:

> 'select()' returned 1:1 mapping between keys and columns Warning
> message: In bitr(x4, fromType = "SYMBOL", toType = "ENTREZID", OrgDb =
> Rat) :
> 98.21% of input gene IDs are fail to map...

Could anyone tell me what happened and how should I do?

Many thanks,
Stanley

Answer: Hi Stanley,

Please try something like this, using the biomaRt package:

require(biomaRt)
mart <- useMart("ENSEMBL_MART_ENSEMBL")
mart <- useDataset("hsapiens_gene_ensembl", mart)
getBM(mart=mart, attributes=c("hgnc_symbol","entrezgene_id"),
filter="hgnc_symbol", values=x, uniqueRows=TRUE)
hgnc_symbol entrezgene_id
1 ANXA4 307
2 CDA 978
3 COL4A1 1282
4 COL4A2 1284
5 COMP 1311
6 CPM 1368
7 CRYAB 1410
8 CRYM 1428
9 CST6 1474
10 DEFB1 1672
11 EEF1A2 1917
12 ERBB2 2064
13 FGFR3 2261
14 FXYD2 486
15 GGT1 2678
16 GGT1 728441
17 GGT1 102724197
18 GLRX 2745
19 GPC3 2719
20 GPX3 2878
21 HCLS1 3059
22 HGD 3081
23 HMOX1 3162
24 HSPA2 3306
25 IGFBP1 3484
26 IGFBP3 3486
27 IL6 3569
28 ITM2B 9445
29 KDELR3 11015
30 LAMB1 3912
31 LAMC1 3915
32 LBP 3929
33 MAL 4118
34 MCAM 4162
35 MITF 4286
36 MYOC 4653
37 NDRG1 10397
38 NID2 22795
39 NNMT 4837
40 NR1H4 9971
41 ORM1 5004
42 PCCA 5095
43 PDXK 8566
44 PDXK 105372824
45 PLS3 5358
46 PTHLH 5744
47 PTPRM 5797
48 PVR 5817
49 RBP4 5950
50 SLC6A12 6539
51 SLPI 6590
52 SOD2 6648
53 SOD2 100129518
54 SPP1 6696
55 STC1 6781
56 TCN2 6948
57 TFPI2 7980
58 TIMP2 7077
59 TOB1 10140
60 WARS 7453

Regarding annotation conversions, there are usually one-to-many mappings, which can complicate things Also important to be aware that many genomic loci are still not well annotated, whilst at others transcripts overlap each other (having different promoters and transcription start sites).

Kevin
Question: Hi guys,

I have a specific problem about using awk or sed to split a big file to different files.
The big file is like this format(3 columns):

C SRR1_45/1 data...
U SRR2_34/2 data...
U SRR1_33/2 data...
C SRR3_22/1 data...
....

I want to extract lines with SRR1 to SRR1.txt, lines with SRR2 to SRR2.txt ...lines with SRRn to SRRn.txt. And the output lines should remove 'SRRi_' symbol. But we don't how many n are there.

e.g. SRR1.txt will contain:
C 45/1 data...
U 33/2 data...

I know it's easy to write a python or perl script to do it. But is there a shell way to do it? taking the advantages of awk or sed.
Let me add some details:
I have 10 such big files to be extracted. And each has more than 1000M lines. So I need to find a efficient way.
The n is random which is not from sequential array.

Thanks!
Tao

Answer: Here is a simple way to do it without sorting and with `awk`:

$ awk '{ split($2, a, "_"); print $1"\t"a[2]"\t"$3 >> a[1]".txt"; }' foo.txt

The file `foo.txt` is a three-column tab-delimited text file containing your data.

Using the `>>` operator appends a line to whatever `SRR*.txt` file exists. Therefore, if you re-run this one-liner, you must first delete any previously-made `SRR*.txt` files, or you will get duplicate lines.

This should be pretty fast, as you're not sorting on IDs. It would be faster, probably, to use a Perl-based approach that opens a pool of file handles, but this should work fine.

Further, if you don't care about the order of lines in the split files, you could use [GNU Parallel][1] with this one-liner to split multiple files `foo1.txt`, `foo2.txt`, etc. simultaneously. Doing the work in parallel may hit a file I/O bottleneck but could give you an overall speed boost, if you use SSDs or other fast storage.

[1]: https://www.biostars.org/p/63816/
Question: I downloaded an SRA file from GEO FTP. I then ran the following command trying to grab the fastq files, but I keep coming up with this error and I can't find anything about it on google.

Here's the command used:

fastq-dump --split-3 '/path/to/SRR1460694.sra'

Error:

fastq-dump.2.5.7 err: connection inconsistent while validating within network system module - failed

I am behind a network firewall, but I have my proxy set up in by .bash_profile so I don't usually ever run into problems downloading, or updating software. I'm not sure what to do at this point. Any ideas?

Answer: If you have to use a HTTP proxy, you have to configure it explicitly for sratools. Sratools do not care about standard unix environment variables like `http_proxy`. NCBI programs have their very own ways of configuration and usage, which often resembles standards from the previous millennium when Fortran programs were fed into computers with punch cards. It seems as if NCBI is still developing with PDP11 in mind.

To make sratools use a HTTP proxy, append the following to your local configuration file which is found at `~/.ncbi/user-settings.mkfg`:

```
/http/proxy/enabled = "true"
/http/proxy/path = "proxy.example.com:8080"
```

See also http://seqanswers.com/forums/showpost.php?p=179416&postcount=8

Support for HTTP proxies has been added to sratools in 2015. Therefore the above settings will not work with some older versions.

It seems that you have already downloaded the SRA file to your local disk. But fastq-dump still needs access the internet in order to download the related reference sequence. Unfortunately, most SRA files are very similar to BAM files, they comprise a mapping of reads to a reference sequence. The reference sequence is not included in the SRA file. SRA files are not self-contained, you need additional external information to read their content.
Question: Hi there, I have a tsv file from GWAS-Catalog like this one in tsv format:
ftp://ftp.ebi.ac.uk/pub/databases/gwas/releases/2018/11/05/gwas-catalog-associations_ontology-annotated.tsv

However, for many annotation tools, they require input data as vcf format.
Therefore, I would like to convert those SNPs into vcf format.
Is there any tools I can use to directly do that ?

Answer: use awk to generate a header and create the VCF field. Here is a basic script.

wget -q -O - "ftp://ftp.ebi.ac.uk/pub/databases/gwas/releases/2018/11/05/gwas-catalog-associations_ontology-annotated.tsv" | awk -F ' ' 'BEGIN{printf("##fileformat=VCFv4.2\n#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n");}/^DATE ADDED/{next}{printf("%s\t%s\t%s\tN\t.\t.\t.\t.\n",$12,$13,$37);}'

##fileformat=VCFv4.2
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
6	32251212	GCST001156	N	.	.	.	.
6	32441753	GCST001156	N	.	.	.	.
6	33075103	GCST001156	N	.	.	.	.
6	32623148	GCST001156	N	.	.	.	.
6	31039078	GCST001181	N	.	.	.	.
6	31125810	GCST001181	N	.	.	.	.
6	31168676	GCST001181	N	.	.	.	.
6	31440051	GCST001181	N	.	.	.	.
Question: Hi,

Does anyone have a strategy for figuring out which restriction enzyme was utilized in a HiC experiment based on the Illumina reads alone?

Thanks

Answer: Figured it out by checking the overrepresented k-mers at the FASTQC report for my data.
Question: After performing gene co-expression network analysis in a case-control study, I am ending up with a bunch of genes whose expression, along with its partners, is dysregulated in two conditions. These highly connected genes are considered as putative "hub" biomarkers. Since, it is not possible for us to consider all of them for experimental validation (for inhibition), what other analysis/tests can be performed to choose one or few of them? I know there is no so caller "master regulator", but somehow if I need to pick only "interesting" genes, what other statistical analysis can be performed?

All suggestions and comments are appreciated.

Answer: Hi vinay.plp,

Since you have a network, you might have a look at topological indexes, like the "betweenes centrality" or the "pairwise disconnectivity index" to rank your genes. As Sam mentioned in his comment, you can also apply the Google page rank algorithm to your network. Have than a look at the top-ranked ones.

You can use external databases to enrich your gene set (e.g. Omim or GO) and filter out classes of interests.

Cheers,

Michael
Question: I have five different cell type i m doing both rld and vsd transformation .I want the vsd data to perform other analyiss such as Self organising maps as well as WGCNA . How can i take out the vsd transformed data

Im using this code

vsd <- varianceStabilizingTransformation(dds, blind=FALSE)

To get all my normalised data set im using this code

resOrdered <- res[order(res$padj),]
summary(resOrdered)
resdata <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized=TRUE)), by="row.names", sort=FALSE)
names(resdata)[1] <- "Gene"

head(resdata)
dim(resdata)
## Write results
write.csv(resdata, file="Primary_cell_line__normalised_data.csv",row.names = FALSE,quote = F)

But this is for the rlog transformation what i saw in the manual if im not wrong , how can i do the same for vsd transformed data ?Any suggestion or help would be highly appreciated

Answer: Your second chunk of code is not for rlog transformed data but for the normalised counts. If you want the data in the same shape as the normalised counts but for rlog or vst transformed data then you need something along the lines of:

vsd <- assay(varianceStabilizingTransformation(dds, blind=FALSE))
rlog <- assay(rlogTransformation(dds, blind=FALSE))
Question: I have a few new Locus tags as follows:

RS09560
RS10020
RS10595

and I need to convert them into old Locus tag which are

RS09560	NA
RS10020	1984
RS10595	2097

NA are the ones with not available old locus tag. I couldn't find any database to do so. And the information available for old locus tag and new locus tag is in the NCBI website

https://www.ncbi.nlm.nih.gov/nuccore/NC_000962.3

I tried to scarp the data from the website using R (rvest) and had written a few lines, but skeptical about the **HTML node** I should extract from the web page **to get the information about old locus tag and new locus tag**.

> library(xml2)
> library(rvest)
> url <- 'https://www.ncbi.nlm.nih.gov/nuccore/NC_000962.3'
> webpage <- read_html(url)
> sb_table <- html_nodes(webpage, ' ')

Please suggest a way to do and also which node could be used to fetch the locus tags.

Answer: I'm a python type of guy, and what I would do is to set up a script to [parse through the genbank file using biopython.][1]

I wouldn't expect to find any source that has any old locus names that are missing from the NCBI genbank files. Gene annotations are constantly updated and it's not unusual for 'novel' genes to pop up, which will not have an old locus name.

Assuming that what you have is a neat list of accessions of interest, it is just a matter of looping through the 'Features' and looking for the accessions, grabbing the old locus names if available, and print a simple list of new vs old accessions.

Even if you have not touched python yet, it is a fairly easy and decent introduction to it if you are willing to give it a go.

[1]: http://biopython.org/DIST/docs/tutorial/Tutorial.html