cannin/biostars_qa huggingface dataset:
Question: Hello,
Cannot find the answer anywhere online. Trying to compile a Perl Module, Math-CDF, for an CLIP (RNA-protein crosslink analysis) and I am getting an error when I run the make commands.

Within the folder I have to run:

$perl Makefile.PL
$make
$make test
$make install

running $perl Makefile.PL runs fine but $make hits the following error:

Skip blib/lib/Math/CDF.pm (unchanged)
cd cdflib && /usr/bin/make LIBPERL_A="libperl.a" LINKTYPE="dynamic" OPTIMIZE="-O3" PREFIX="/usr/local/ActivePerl-5.24/site" PASTHRU_DEFINE=' ' PASTHRU_INC=' '
make[1]: Nothing to be done for `all'.
Running Mkbootstrap for CDF ()
chmod 644 "CDF.bs"
"/usr/local/ActivePerl-5.24/bin/perl-dynamic" -MExtUtils::Command::MM -e 'cp_nonempty' -- CDF.bs blib/arch/auto/Math/CDF/CDF.bs 644
/usr/bin/gcc -c -fno-common -DPERL_DARWIN -no-cpp-precomp -mmacosx-version-min=10.9 -arch x86_64 -DUSE_SITECUSTOMIZE -DPERL_RELOCATABLE_INCPUSH -fno-strict-aliasing -pipe -fstack-protector -DPERL_USE_SAFE_PUTENV -O3 -DVERSION=\"0.1\" -DXS_VERSION=\"0.1\" "-I/usr/local/ActivePerl-5.24/lib/CORE" CDF.c
gcc: error: unrecognized command line option '-no-cpp-precomp'
make: *** [CDF.o] Error 1

Any ideas? I've tried removing the -no-cpp-precomp flag but still doesn't run.
P.S I am running a Mac O.S Yosemite Version 10.10.5 and I working with Perl modules for the first time.
I have make and gcc under:
/usr/bin/gcc
/usr/bin/make
Thanks.

Answer: Can't seem to delete the post, but thanks everyone for your help. The solution here was just to delete the -'-no-cpp-precomp' flag from the makefile and re run. It is running smoothly now, I am using it for this bioinformatics protocol:
https://zhanglab.c2b2.columbia.edu/index.php/CTK_Documentation
Question: I am using biomart to convert multiple gene names to Entrez IDs, however, I found that for some gene names, it returned multiple IDs. Shouldn't the name and ID have a one to one correspondence? I also tried to use all the returned IDs and search in NCBI, but it seems only one of them is the same with my inquiry. For example, try "POTED" in biomart, it returns three IDs - "388468,339010,317754". But only searching 317754 in NCBI will give you the right gene symbol POTED. Why biomart returns multiple IDs?

Answer: Gene names aren't unique. The closest to an actually unique ID that I normally use would be Ensembl gene IDs (this is why I normally recommend using Ensembl annotations, when available).

Edit: I'll add that a single Entrez ID can have multiple gene names (on more than one chromosome or even on different strands of the same chromosome). Yeah, that can make life difficult.
Question: Does anyone know of a database of integration sites of known Human Endogenous Retroviruses (HERVs). I've come across [this resource][1] which came from [this paper][2] in 2003 but the data behind the website is missing.

I could probably take their method and scan the genome myself, but a citable database would make everything easier.

[1]: http://herv.img.cas.cz/
[2]: http://nar.oxfordjournals.org/content/32/suppl_1/D50.abstract

Answer: Your best bet might be the rmsk table in UCSC

Repeats -> Repeatmasker -> rmsk

http://genome.ucsc.edu/cgi-bin/hgTables

the type column will show "ERV"
Question: Would you advice how to get promoter sequences for all or many human genes - in flat file(s) or by SQL query? I understand there can be multiple definitions for promoter region, but just anything universal would work.

Answer: My best bet would be to use BioMart's Martview, select a database, filter by the gene IDs you have (there are other ID options there too), and the use the sequence option in the attributes to determine which parts of the gene you want, be exon, intron, promoter, upstream, downstream, etc.

I used this tool to get many upstream regions for mouse genes using just NCBI's gene IDs as input.
Question: I am trying to run secretomep for which psort II is listed as a dependancy package. However, I cannot find any copies of a standalone version of psort II for linux ubuntu anywhere.. Not one that is free anyway. I found a company ( [http://www.genscript.com/psort.html][1] ) offering the tool for a price. Since, I am doing this work independantly and have no funding, I can't really afford to get commercial tools. But I really need secretomep to work too. Can anyone help me out with a copy of PSORT II and WoLFPSORT if you already have it, or please direct me to a page where I could download it.

[1]: http://www.genscript.com/psort.html

Answer: Well it looks like someone has put up code on GitHub for Wolf-pSort as the original webserver is now defunct [Wolf-PSort Github][1]. PSortII looks like it may no longer be easily available for local install. The email address for the developer is knakai@ims.u-tokyo.ac.jp but I don't know if they are still active (may have retired), etc. Hopefully someone may have a locally downloaded copy of the source code that they can share. If someone DOES have it, and if the license allows (it should, I believe it was originally under the GPL, I would suggest creating a GitHub repository like someone did for Wolf-PSort to keep it available and distributed.

[1]: https://github.com/fmaguire/WoLFPSort
Question: Hi all,

I have expression measure of genes (two sets: target and contro genes) as measured by RNAseq, which is computed in tabular matrix.

head(A)
V1 V2 V3 V4 V5 V6
1 chr20:38788661-38812050_ENSDART00000136771 1.912 1.912 1.912 4.770 0.703
2 chr4:13971999-13980660_ENSDART00000067040 0.502 0.502 0.502 1.242 0.011

V7 V8 V9 V10 V11 V12 V13
1 0.137 6.534 4.689 7.086 2.646 4.376 0.684
2 0.007 0.684 1.125 2.819 2.392 1.406 0.364
head(B)
V1 V2 V3 V4 V5 V6
1 chr19:24974067-24977991_ENSDART00000052421 0.005 0.005 0.005 0.025 0.139
2 chr8:26148102-26150871_ENSDART00000049793 2.312 2.312 2.312 6.153 1.109
3 chr8:47908733-47913950_ENSDART00000025620 2.446 2.446 2.446 9.660 5.411
V7 V8 V9 V10 V11 V12 V13
1 0.807 0.010 0.232 5.963 0.022 0.149 3.162
2 0.291 4.556 3.214 16.486 2.158 5.599 0.939
3 7.745 11.614 13.150 51.680 8.557 18.356 36.378

PLOT:

par (mfrow=c(1,2))
plot (log2(A$V2/A$V5), log2(A$V8/A$V11), xlab="log2(FC)", ylab="log2(FC)", xlim=c(-5,5), ylim=c(-5,5), col="red")
abline(v=0, lty=2); abline (h=0, lty=2)
plot (log2(B$V2/B$V5), log2(B$V8/B$V11), xlab="log2(FC)", ylab="log2(FC)", col="blue", xlim=c(-5,5), ylim=c(-5,5))
abline(v=0, lty=2); abline (h=0, lty=2)

How could i plot both in the same plot?

cheers

Chirag

Answer: a)

plot(c(log2(A$V2/A$V5),log2(B$V2/B$V5)), c(...), col=c(rep("red", nrow(A)),rep("blue", nrow(B)))

b)

with ggplot2 its much more elegant. You would have to create one data frame with a structure like:

gene name x y type
gene1 3.2 1.2 control
....
geneN 1.2 3.2 target

and then `ggplot(df, aes(x=x,y=y, colour=type))+geom_point()`
Question: Is there a way to add attributes to the SVG output of the ETE tree drawing toolkit?

I am drawing phylogenetic trees as svg and would like to add the NCBI taxonomy ID as ID attribute to the node in the SVG. This would let me select them easily with java script, when using the SVG in a html page.

Answer: Unfortunately this is not possible. The closest you can get is to generate a png image and use the pixel map dictionary returned by the tree.render() function.
Question: When I input the word *titin* in [UCSC genome][1] I get multiple relevant hits as shown [here][2]. My interest in defining the DNA coding region of the titin gene. Specifically when I click on one of my search hits

[TTN (uc021vtb.4) at chr2:178525989-178807423 - Homo sapiens titin (TTN), transcript variant N2-B, mRNA. (from RefSeq NM_003319)][3]

and view the DNA, it is positioned at **chr2:178,525,989-178,807,423**. My guess is that this is the DNA coding strand.

Now my confusion comes when I refer to the same titin gene from [NCBI RefSeq][4]. The sequence is **chr2:178525989..178807423, complement**. And when I look at the exact NCBI's FASTA sequence, it is the reverse complement of that from UCSC Genome. So which is the correct coding DNA?

[1]: https://genome.ucsc.edu/cgi-bin/hgGateway
[2]: https://genome.ucsc.edu/cgi-bin/hgTracks?hgsid=647507939_GytALTLl6nFjgkaznatQhwHF8hTg&org=Human&db=hg38&position=titin&pix=921
[3]: https://genome.ucsc.edu/cgi-bin/hgTracks?db=hg38&lastVirtModeType=default&lastVirtModeExtraState=&virtModeType=default&virtMode=0&nonVirtPosition=&position=chr2%3A178525989-178807423&hgsid=647507939_GytALTLl6nFjgkaznatQhwHF8hTg
[4]: https://www.ncbi.nlm.nih.gov/gene/7273

Answer: We always store our coordinates in relation to the positive strand, because it makes coordinate math easier to deal with. When you click on the transcript in question, you can see the following information:
Position: hg38 chr2:178,525,989-178,807,423 Size: 281,435 Total Exon Count: 191 Strand: -

The position here indicates only the genomic range of the item, and has nothing to do with the strandedness of the transcript.

When you choose the DNA output link from the details page for your transcript, the DNA returned is reverse complemented. You can verify this by taking the first set of capitalized bases (the first exon) and BLATting the sequence against hg38, where you can see that the first exon is at the "rightmost" end of the transcript:

http://genome.ucsc.edu/cgi-bin/hgTracks?hgS_doOtherUser=submit&hgS_otherUserName=chmalee&hgS_otherUserSessionName=hg38_ttn_first_exon

For more information on our coordinate system, please see the following blog post:

http://genome.ucsc.edu/blog/the-ucsc-genome-browser-coordinate-counting-systems/

Lastly, when you share links to the Genome Browser do not include the "hgsid" in the URL, as the hgsid indicates session specific information to our servers, and thus anyone who uses that hgsid and browses to a new location, or makes custom tracks, etc will cause the link you sent to not be what you intend. If you want to create a shareable link, please see our Session documentation:

http://genome.ucsc.edu/goldenPath/help/hgSessionHelp.html

If you have further questions about UCSC data or tools feel free to send your question to one of the below mailing lists:

- General questions: genome@soe.ucsc.edu
- Questions involving private data: genome-www@soe.ucsc.edu
- Questions involving mirror sites: genome-mirror@ose.ucsc.edu

ChrisL from the UCSC Genome Browser
Question: Dear All,

I am doing few histone protein's SNP analysis.

I have used SIFT and Polyphen for my initial analysis. Now I would like predict its functional impact.

Can anyone please suggest me few tools or software for the prediction of the functional

effect of nsSNPs. Thank you

Answer: You can try [SnpEff][1] + SnpFit (dbNFSP)

NOTE: RamRs noted that the input data format is NOT from a sequencing experiment, but another format that lacks of the corresponding genomic coordinates. In such case, this recommendation is not appropriate for the analysis.

[1]: http://snpeff.sourceforge.net/
Question: Dear members,

I am using "An end to end workflow for differential gene expression using Affymetrix microarrays" as described [here][1]
at one of the step it requires to map probe ids to gene names and symbols.

The methods uses

anno_palmieri <- AnnotationDbi::select(hugene10sttranscriptcluster.db,
keys = (featureNames(palmieri_manfiltered)),
columns = c("SYMBOL", "GENENAME"),
keytype = "PROBEID")
However, I get an error, saying:-

Error in .testForValidKeys(x, keys, keytype, fks) :
None of the keys entered are valid keys for 'PROBEID'. Please use the keys method to see a listing of valid arguments.

On some digging I found that:-

my keys look like this:

head(featureNames(palmieri_manfiltered))

> "1007_s_at" "1053_at" "117_at" "121_at" "1255_g_at" "1294_at"

whereas if I look at the key type probe id from hugene10sttranscriptcluster.db using

head(keys(hugene10sttranscriptcluster.db, keytype = "PROBEID"))

it looks like this

> [1] "7892501" "7892502" "7892503" "7892504" "7892505" "7892506"

Can anyone help me in fixing this issue or put me in right direction. Thank you

[1]: https://www.bioconductor.org/packages/devel/workflows/vignettes/maEndToEnd/inst/doc/MA-Workflow.html#11_annotation_of_the_transcript_clusters

Answer: Which array are you using? Your probe IDs are not from the Human Gene 1.0 ST Array - they seem to be from the HG-U95. So, you could try either of these DBs instead:

hgu95av2.db
hgu95a.db

Kevin
Question: I have different value for test$pval and and other.pval that should return the same log-rank test p-value, no ?

surv_object <- Surv(time = as.numeric(final_group_annotated[,time]), event = final_group_annotated[,opt$endPoint] )it1 <- survfit(surv_object ~ group, data = final_group_annotated)

test=surv_pvalue(fit1, final_group_annotated)
print(test$pval)

> 0.3427874

fit.coxph <- coxph(surv_object ~ group, data = final_group_annotated)
print(summary(fit.coxph))

coxph(formula = surv_object ~ group, data = final_group_annotated)

n= 188, number of events= 30

coef exp(coef) se(coef) z Pr(>|z|)
groupgroup2 -0.3482 0.7060 0.3689 -0.944 0.345

exp(coef) exp(-coef) lower .95 upper .95
groupgroup2 0.706 1.417 0.3426 1.455

Concordance= 0.538 (se = 0.049 )
Likelihood ratio test= 0.87 on 1 df, p=0.4
Wald test = 0.89 on 1 df, p=0.3
Score (logrank) test = 0.9 on 1 df, p=0.3

other.pval <- coef(summary(fit.coxph))[,5]
print(other.pval)

> 0.3452066

Answer: The p-value returned by `Surv()` is the log-rank p-value from the score test; so, you should be comparing this to the log-rank p-value from the Cox PH model. I went over this recently in this thread: https://www.biostars.org/p/401657/#401759

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

The way that you are doing it, i.e.:

coef(summary(fit.coxph))[,5]

This will return multiple p-values when there are multiple strata, one p-value representing a single stratum / coefficient compared to the reference stratum.

To extract the log-rank p-value for the overall model, you need to use:

summary(coxfit)$sctest[3]

Kevin
Question: Hi,

I have two data frames, both containing the columns `GeneName` and `logFC`.

Using `intersect()`, I have obtained the common genes between the two. I want to know how to merge my two dataframes based on the common dataframe.

Thank you.

Answer: **1. Get common column names**

common_column_names <- intersect(names(dataframe_1), names(dataframe_2))

**2. Use this character vector as your `by` parameters in the `merge` function.**

merge(dataframe_1, dataframe_2, by=common_column_names, all.dataframe_1=TRUE)

** P.S - Your title of the question should be "How to merge 2 dataframes based on common column names?" I have edited the title.
Question: I have thousands of EST sequences in a fasta file. How to subset a set of sequences based on sequence ID or name using R?

Sequence example:

>gi|296783888|gb|GW992815.1|GW992815 UAS-Mi10 Complementary DNA of mulberry (Morus indica) Morus indica cDNA 5' similar to Putative phosphoribosyltransferase/phosphoribosylanthranilate-like gene from Morus indica, mRNA sequence

GCAGCCGTCGGATCGTGAGCGTGATCGCGTGGCTAGTCGGGTTGGCGAAATGGTTGGATGATATCCGGAG

GTGGAGGAACCCCATTACCACGGTATTGGTCCACATCTTATATTTAGTGCTTGTTTGGTACCCGGATTTG

ATTGTCCCAACCGGGTTTTTATATGTGTTCCTAATCGGTGTATGGTACTATCGGTTTCGGCCCAAGATAC

CAGCGGGTATGGATACCCGACTCTCACAAGCTGAAGCGGTTGACCCGGATGAGCTTGATGAGGAATTCGA

CACCATACCGAGCTCAAAACCACCCGACATAATCAGGGTCCGGTATGACCGGTTGCGGATATTGGCAGCC

CGGGTTCAAACGGTTTTGGGTGATTTTGCAACACAAGGGGAGCGGGTTCAGGCCTTGGTTAGCTGGAGGG

ACCCAAGGGCCACAAAATTGTTCATAGGCGTGTGCTTGGCCATAACAATAATTCTCTATGTGGTGCCACC

CAAAATGGTTGCCGTGGCACTTGGATTCTACTATTTACGACACCCCATGTTCCGAGACCCCATGCCTCCT

GCAAGCTTGAATTTCTTCAGAAGGCTTCCAAGCCTTTCAGACCGCTTTAATGTAGATTAGAATATTATAT

GATTATTAGTAGGCCCAA

>gi|296783887|gb|GW992814.1|GW992814 UAS-Mi9 Complementary DNA of mulberry (Morus indica) Morus indica cDNA 5' similar to Dehydration-responsive protein RD22, Similar to BURP domain-containing protein like gene from Morus indica, mRNA sequence

AAGCAGTGGTCTAGAACCAGAGTGGCCCCTGCGATGCAGGTATCATCTCTATTATCAAAAGGGATAAGGG

GTGGATCCGTCGGGGATTTGAGTCTCACATGGTCGCTGATAACTTATTGAATGGATATTGGATTGTGTGC

AGTGCGACCTAAACAGGATTGCCGTTGGGGCCTGTGGTCAGAGATACCCCACACTTCTCAACTCCCAAAT

TGGATCTTGTTCCTTGTTTTCCTGTATTAAGCCTGACCCCTGAGGCTTTCGCCACTGCCAACTGGGTGCC

GCCTGCTGACTTCTGATTCCCCGTGCTAACGGTTACTCCCGATTCCTTATCCACATCGAAGATGAACTAT

TGACTTCCGCAAACTCAAAAGGCTGCAAGATATCACTGACCGCTGTCGGGATCCGCGATCGGCATATACG

CGAAATCCGATCCCGGATCCCGGGACTGCAGACGGCTGAA

Like using this header:

>gi|296783888|gb|GW992815.1|GW992815 UAS-Mi10 Complementary DNA of mulberry (Morus indica) Morus indica cDNA 5' similar to Putative phosphoribosyltransferase/phosphoribosylanthranilate-like gene from Morus indica, mRNA sequence

or by using gi number?

How to do this in R?

Answer: library(Biostrings)
f <- readDNAStringSet("sequences.fa")

You can then use the names(f) accessor to either directly match the names of the sequences, grep with gi accession, or even just split the names and directly match the gi accessors. The f object can be directly subset (e.g., f[c(1,3,5,6)]).
Question: Does anyone know of a reliable package for doing pairwise DE analysis on microarray data in the form of an ExpressionSet object? I've used edgeR before for doing pairwise comparison with rna-seq, but I'm a little new to this and can't really find anything to do it on affy data specifically...

Answer: I use *samr*:

R. Tibshirani, G. Chu, Balasubramanian Narasimhan and Jun Li (2011). samr: SAM: Significance Analysis of Microarrays. R package version 2.0. http://CRAN.R-project.org/package=samr

Expression data are a matrix E, with identifiers PROBE_IDENTIFIERS and symbols GENE_SYMBOLS, and labels indicating pairs. Labels are integers starting at 1, with positive or negative indicating paired samples. Works for any array (or indeed any matrix of numbers).

```
data = list(x=data.matrix(E),
y=labels,
geneid=PROBE_IDENTIFIERS
genenames = GENE_SYMBOLS,
logged2=TRUE)
samr.obj = samr(data, resp.type="Two class paired", nperms=100)
```
Question: Hi,

I received some bams where duplicates were marked using Picard. I'd like to unmark all duplicate reads for a study we're doing. Does anyone know of an existing tool? I haven't found anything in google land.

If there isn't an existing tool, could you provide some guidance how to unmark them myself? I've read [the SAM spec][1] and understand I need to change the bit flag, but I'm not sure the best way to approach it.

Thanks!

Mark

[1]: http://samtools.github.io/hts-specs/SAMv1.pdf

Answer: as said Devon, writing this tool is very easy ;-)

I wrote it and put it here: https://github.com/lindenb/jvarkit/wiki/Biostar106668

Example:

```
$ curl -s "ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/phase1/technical/ncbi_varpipe_data/alignment/NA12751/NA12751.chrom20.ILLUMINA.mosaik.CEU.low_coverage.20101123.bam" |\
samtools view -h - |\
head -n 10000 |\
java -jar dist/biostar106668.jar -b > out.bam
```
Question: I like to know is there any way to merge the same coordinate in a table.

Table:

chr1	155944562	155945214	fantom_neuron	GSM1554667	9.84447
chr1	155944562	155945214	fantom_neuron	GSM1554672 7.43630
chr1	155944562	155945214	fantom_neuron	GSM1554678	32.77627
chr1	155945743	155946196	fantom_neuron	BAMPE	3.87072
chr1	155945743	155946196	fantom_neuron	GSM1554666	18.14939
chr1	155945746	155946202	fantom_neuron	GSM1554655	1.14939

Expected table:

chr1	155944562	155945214	fantom_neuron	GSM1554667,GSM1554672,GSM1554678	9.84447,7.43630,32.77627
chr1	155945743	155946196	fantom_neuron	BAMPE,GSM1554666	18.14939 ,3.87072
chr1	155945746	155946202	fantom_neuron	GSM1554655	1.14939

Answer: Agreeing with Wouter, check out [bedtools merge](https://bedtools.readthedocs.io/en/latest/content/tools/merge.html):

bedtools merge -i foo.bed -c 4,5 -o collapse
Question: I have the following information from blastx annotations of bacterial genes predicted by prodigal:

Sequence name	Sequence desc.	Sequence length	Hit desc.	Hit ACC
gene_1_contig_1	excinuclease ABC subunit A	228	gi|1055624747|ref|WP_067265422.1|excinuclease ABC subunit A [Sulfitobacter sp. HI0054] gi|1024544140|gb|KZY51396.1| excinuclease ABC subunit A [Sulfitobacter sp. HI0054]	WP_067265422, KZY51396
gene_2_contig_1	excinuclease ABC subunit A	210	gi|1055651942|ref|WP_067291557.1|excinuclease ABC subunit A [Sulfitobacter sp. EhC04] gi|1032103716|gb|OAN76192.1| excinuclease ABC subunit A [Sulfitobacter sp. EhC04]	WP_067291557, OAN76192
gene_3_contig_1	MFS transporter	432	gi|1055624744|ref|WP_067265419.1|MFS transporter [Sulfitobacter sp. HI0054] gi|1024544139|gb|KZY51395.1| hypothetical protein A3734_05250 [Sulfitobacter sp. HI0054]	WP_067265419, KZY51395
gene_4_contig_1	MFS transporter	561	gi|1055624744|ref|WP_067265419.1|MFS transporter [Sulfitobacter sp. HI0054] gi|1024544139|gb|KZY51395.1| hypothetical protein A3734_05250 [Sulfitobacter sp. HI0054]	WP_067265419, KZY51395

I wish to fetch gene symbols using the information (either the gi identifiers or the protein accessions) from the blastx results; may be using entrex efetch.

**So, the result would be as below:**

Gene Name Gene symbol
excinuclease ABC subunit A UvrA

See, the link [here][1]. However, I am not sure how to proceed in this case. Can anybody please suggest something?

[1]: https://www.ncbi.nlm.nih.gov/gene/?term=UvrA%20-Pppraxal%20%3Ap[rajapati026%40gmail.com

Answer: Hi,

You can use GI ids to retrieve associated information from uniprot "Retrieve/ID mapping" [UniProtKB][1].
Here, from "GI number" to "UniProtKB" should be selected and it will give output with all the information you need in tabular format and you can select columns of your interest.

Another way is to use [batchentrez][2] to get gene bank data and you need to parse the information.

[1]: http://www.uniprot.org/uploadlists/
[2]: https://www.ncbi.nlm.nih.gov/sites/batchentrez
Question: I was going to compare gene expression in primary tumor and normal control.

So I downloaded PAN cancer data from xena.USCS.edu (https://tcga.xenahubs.net/download/TCGA.PANCAN.sampleMap/PANCAN_clinicalMatrix
for clinical matrix and
https://toil.xenahubs.net/download/tcga_RSEM_gene_fpkm
for gene expression TOIL Gene RSEM fpkm dataset.

Then I selected patient that have matched normal and primary tumor sample types, and I found that for each gene, fpm from normal and primary tumor are JUST THE SAME!

I wonder what's wrong with this.

Is there anybody helping me? Thanks in advance!

Answer: In almost all cases cases, each TCGA patient only had tumor RNAseq done. ([see previous discussions][1]). Are you sure that there was actually normal RNAseq available, or have you pulled down the same file twice through different linkages?

[1]: https://www.biostars.org/p/179681/#179746
Question: I've some snps covering my gene X. I want to know if any of those snps is included in a domain of the protein encoded by my gene X. Anyone can suggest me how can I know the coordinates of the domains of my gene to see if my snps are included in these regions? Or any other method to know this?

Thanks in advance,

Answer: [VEP](http://www.ensembl.org/info/docs/tools/vep/script/vep_options.html) has got an option `--domains`

> `--domains` Adds names of overlapping protein domains to output. Not used by default
Question: Hello,

Can anyone suggest a simple R package to perform GO term analysis. Specially I am looking for a package that can produces publication quality figures and heatmaps regarding GO term results.

Thanks in advance

Answer: There are quite a few packages out there, and there was a thread on this already: https://www.biostars.org/p/222081/

See also my answer to this thread on GO term reliability: https://www.biostars.org/p/272473/#272486

One that is not mentioned in these threads is topGO

Another one with heatmaps is GOexpress