cannin/biostars_qa huggingface dataset:
Question: Hello everyone,

I am currently working on reverse engineering Gene Regulatory Networks from RNA-seq data. I have data from different stages of heart development. My question is if you have different transcripts from the same gene in the RNA-seq data (which is ofcourse unavoidable), should one use information from both transcripts to build the network? because obviously the network is at the gene level and therefore using transcript-level information can be tricky because the transcripts will have different expression values. How will this affect the overall network? plus won't the network inference programs get confused as there will be 2 entries with the same gene name but different values? One will have to make the gene names unique and so these will become 2 different nodes in the network. In such a case, is it OK to use this kind of information as people do build networks from RNA-seq data? Should one keep 1 transcript per gene? but how to select one transcript and on what basis?

Any thoughts?

Thanks!!

Answer: First, with RNA Seq data, you can always choose to use the Gene level information or the Transcript level information. Personally I prefer using the gene level information (counting the coverage per gene using tools such as HTSeq) as they are much easier to interpret and with less noise.

Now with the network analysis. I have only used WGCNA for de novo network analysis and it relies heavily on the correlation matrix. So it can either be gene - gene correlation matrix or transcript - transcript correlation matrix. But if you provide the transcript information, you are building the transcription co-expression network which is slightly different.

So, basically, you can give each transcript an unique identifier and perform the network analysis as if it is gene co-expression network analysis or you can just try to obtain the gene level information and perform the gene base analysis.

One important thing to note is that you should keep a very good record as to the data pre-processing procedures that you have taken. As these network co-expression analysis usually relies on clustering, a small change of input e.g. different form of normalization/ different number of genes (transcript) can alter the final result. So if you are choosing to represent a gene by one particular transcript (which I think is very ambiguous and very likely one will not find a gold standard), you should make sure you have a clear rule for replication.

[Here][1] is the WGCNA page

[1]: http://labs.genetics.ucla.edu/horvath/CoexpressionNetwork/Rpackages/WGCNA/
Question: Hi,

I am working on small RNA-Seq data analysis, I have query regarding trimmomatic tool which widely used for adapter trimming from reads.

In the command below:

trimmomatic SE -phred33 input.fastq Out_trimmomatic.fastq ILLUMINACLIP:/trimmomatic-0.38-1/adapters/TruSeq3-SE.fa:2:30:10 LEADING:30 TRAILING:30 MINLEN:36

TruSeq3-SE.fa its given in software package, Can we use same for all RNA-seq data-sets?
What is the purpose of 2:30:10 after the adapter file?

Thanks in Advance..!!

Answer: From the [Trimmomatic Website][1]

> Step options:
>
> ILLUMINACLIP:<fastaWithAdaptersEtc>:<seed mismatches>:<palindrome clip
> threshold>:<simple clip threshold>

> fastaWithAdaptersEtc: specifies the
> path to a fasta file containing all the adapters, PCR sequences etc.
> The naming of the various sequences within this file determines how
> they are used. See below.

> seedMismatches: specifies the maximum
> mismatch count which will still allow a full match to be performed

> palindromeClipThreshold: specifies how accurate the match between the
> two 'adapter ligated' reads must be for PE palindrome read alignment.

> simpleClipThreshold: specifies how accurate the match between any
> adapter etc. sequence must be against a read.

So, you could just use it as is for any [RNAseq] data set with TruSeq2 adapters, but you might want to check your results and then tweak the parameters if it over- or undertrims.

[1]: http://www.usadellab.org/cms/?page=trimmomatic

Edit: reformatted quotation from Trimmomatic website to make it more legible
Edit2: addressed TruSeq
Question: Is there anyway to get BioGrid gene interaction network automatically using R?

Answer: yes, use the REST service http://wiki.thebiogrid.org/doku.php/biogridrest and the R methods to access the web http://www.math.ucla.edu/~anderson/rw1001/library/base/html/read.table.url.html
Question: How do I convert GFF file to a GTF file? Is there any tool available?

Answer: The easiest way is to use the gffread program that comes with the Cufflinks software suite (Tuxedo)

gffread my.gff3 -T -o my.gtf

See `gffread -h` for more information
Question: I have a list of enriched GO terms and would like to have something like the following graph. This graph was done using RAMIGO R package, but now is giving me errors.

Thanks

Here is the RAMIGO error:

```r
> library(RamiGO)

**Loading required package: gsubfn**
**Loading required package: proto**

> d = read.table("same_GOs.txt", header = TRUE) # note that required a header named 'GOs'
> v <- as.vector(d$GOs)
> color <- c("lightblue")
> #v <- as.vector(d$GOs)
> goIDs <- v
> pngRes <- getAmigoTree(goIDs=goIDs, color=color, filename="example", picType="png", saveResult=TRUE)

Error in readPNG(aa) :
unable to open

301 Moved Permanently

Moved Permanently

The document has moved here.

> sessionInfo()

R version 3.0.3 (2014-03-06)
Platform: x86_64-apple-darwin10.8.0 (64-bit)
locale:
[1] C/UTF-8/C/C/C/C

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

other attached packages:
[1] RamiGO_1.8.1 gsubfn_0.6-6 proto_0.3-10

loaded via a namespace (and not attached):
[1] BiocGenerics_0.8.0 RCurl_1.95-4.5 RCytoscape_1.12.0 XML_3.95-0.2
[5] XMLRPC_0.3-0 bitops_1.0-6 graph_1.40.1 igraph_0.7.1
[9] parallel_3.0.3 png_0.1-7 stats4_3.0.3 tcltk_3.0.3

[13] tools_3.0.3
```

Answer: This package seems to try to use the Amigo visualization web service at http://amigo.geneontology.org/visualize to draw the graph but the endpoint was moved to **http://amigo.geneontology.org/amigo/visualize**. You need to update R and the BioC package to the latest version, the version you have is very old, and then try again.
Question: Hi everyone,

I have used the "findOverlaps" function in R to find which positions of my two datasets overlap. I have also used "countOverlaps" to see how many overlaps I have.

What I want to do now, is to find the exact coordinates which overlapped. I was looking how to do it using findOverlaps but it has no such options and by doing a google search, I couldn't find much help.

A little help on this would be greatly appreciated.

Thank you

Answer: The `findOverlaps()` method is not the right tool for the job. Instead, you want the set intersection which can be obtained with `intersect()` on two GRanges objects.

# create two GRanges objects for testing
gr1 = GRanges(seqnames='chr1',ranges=IRanges(start=seq(1,1001,100),width=50))
gr2 = GRanges(seqnames='chr1',ranges=IRanges(start=seq(25,1025,100),width=50))

# and "intersect" them, finding the set of regions
# formed by the overlaps of gr1 and gr2
gr3 = intersect(gr1,gr2)

# and the result
gr3

GRanges object with 11 ranges and 0 metadata columns:
seqnames ranges strand

[1] chr1 [ 25, 50] *
[2] chr1 [ 125, 150] *
[3] chr1 [ 225, 250] *
[4] chr1 [ 325, 350] *
[5] chr1 [ 425, 450] *
[6] chr1 [ 525, 550] *
[7] chr1 [ 625, 650] *
[8] chr1 [ 725, 750] *
[9] chr1 [ 825, 850] *
[10] chr1 [ 925, 950] *
[11] chr1 [1025, 1050] *
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
Question: Hi,

I have a huge list of gene ids which are 0610005C13Rik, 0610007P14Rik etc in mouse RNA-Seq data. I want to convert it into either Refseq id or common gene symbol in R using Bioconductor package. I'm new to this and open to trying new ways.

Answer: You can use [NCBI Datasets](www.ncbi.nlm.nih.gov/datasets/) for this. Go to the [NCBI Data Tables](https://www.ncbi.nlm.nih.gov/datasets/tables/genes/) page and upload a file containing the list of names of the genes. Choose 'mouse' as the organism. The output table can be modified to include/exclude specific columns of interest to you and downloaded as a tab-delimited file.

Alternatively, you can use the [command line tool](https://www.ncbi.nlm.nih.gov/datasets/docs/command-line-start/) to obtain the data in json format and parse it using a tool like `jq` to extract fields of interest.
Question: Does any one know how to extract SRX ids from SRR ids. Thanks.

ex:

***Input***

SRR1294518

***Output***

SRR1294518 SRX549132

Answer: Have a look at the SRAdb package in Bioconductor. You do have to download the SRA database backup (I suspect there's a way to remotely query it, but perhaps not). Then you can use the sraConvert() command. This also allows you to download datasets associated with a run/sample/whatever easily.
Question: Hi,

I am trying to figure out whether my gene of interest is expressed differentially between two tissues. I have been looking at GEO, but I am not clear on how to actually find my gene. I have tried all possible names I can come up with, but GEO does not seem to know the gene.

I guess this leaves me with the following questions:

1. Can I search GEO based on chromosome location?
2. Do you know of other databases that would allow me to do this?
3. Should I download the raw reads and redo the analysis?
4. other ideas?

As always, thanks for any help!

Answer: Are you looking at a particular dataset in GEO? More specific details would be helpful.

You might want to have a look at GTEX, which is a large systematic dataset of gene expression and eQTLs across many human tissues: http://www.gtexportal.org/home/
Question: I am currently using BWA-MEM to map metagenomic reads to a small (~12kb) virus reference genome for a large number of samples. I have a lot of reads per sample - upwards of 10 million.

The majority of the reads (~99%) do not match my reference - I am not interested in analysing these reads. However, the .sam file that BWA-MEM produces stores the unmapped reads, and I am stuck with these huge files that are taking up disk space (I do convert these into sorted .bam files and delete the .sam file, but they're still huge).

Is there any way to stop BWA-MEM from storing the unmapped reads in the .sam file and only keep the mapped reads?

Thank you for letting me know.

Answer: bwa mem (...) | samtools view -F 4 -o onlyMapped.bam

See http://broadinstitute.github.io/picard/explain-flags.html
Question: Hi all!

I am using bwa 7.12 to map reads to reference genome. I was analysing exome when after 40 min of running the error appeared:

[gzclose] buffer error

I've red that I should upgrade the zlib package for actual bwa version, but I have no clue how to do that (to make it compatible with bwa).

In addition, I didn’t have this problem before upgrading my Ubuntu to 16.4 version. Could this upgrade mess with my bwa configurations?

Any help would be much appreciated.

To be more specific I've tried:

$ wget http://downloads.sourceforge.net/project/bio-bwa/bwa-0.7.12.tar.bz2
$ tar xvjf bwa-0.7.12.tar.bz2
$ cd bwa-0.7.12
$ sed -e's#INCLUDES=#INCLUDES=-Izlib-1.2.11/ #' -e's#-lz#zlib-1.2.11/libz.a#' Makefile > Makefile.zlib
$ wget http://zlib.net/zlib-1.2.11.tar.gz
$ tar xvzf zlib-1.2.11.tar.gz
$ cd zlib-1.2.11
$ ./configure
$ make
$ cd ..
$ make -f Makefile.zlib

but it didn't work...

Best,

Agata

Answer: First, I would suggest you check your fastq files are not truncated or corrupted.

If you bwa is from Ubuntu repositories, then the error is probably not related to an old zlib, as Ubuntu Xenial zlib is 1.2.8, and the error you reported is linked to zlib older than 1.2.4. Anyway, you may install a newer bwa:

wget https://sourceforge.net/projects/bio-bwa/files/bwa-0.7.15.tar.bz2
tar -xvjf bwa-0.7.15.tar.bz2
cd bwa-0.7.15/
make

Then add this folder to the PATH.
Question: What are some good resources to learn shell script for NGS pipeline development?

How much shell script should one know to develop an intermediate level pipeline for NGS data analysis?

Can someone suggest some good resources, tutorials?

Answer: I use http://www.ruffus.org.uk/ to develop python pipelines for NGS data and it is a great library.
Question: Hi,

I was wondering how databases as 1000genomes and ESP6500 deal with allele frequencies for positions on the X chromosome. Do they only take females to calculate the number of alleles & allele frequencies or do they take into account that a male only has one X?

Answer: They just count how many of each allele exist overall.

Let's say we have a variant on X and a group of two men, two women. Their genotypes are A, G, AG and GG - that's six alleles, two A and four G. We would call A as the minor allele with a frequency of 0.33.

[Here][1] is an example. Look at the top line of the table, ie ALL. There are 384 C alleles and 1275 T alleles in the total population, so the MAF is 0.231 (384/1659 = 0.231). The genotypes are 129 C men, 396 T men, 37 CC women, 181 CT women and 349 TT women. C frequency = 129 (C men) + 2 x 37 (CC women) + 181 (CT women) = 384. T frequency = 396 (T men) + 181 (CT women) + 2 x 349 (TT women) = 1275.

[1]: http://www.ensembl.org/Homo_sapiens/Variation/Population?db=core;r=X:84436603-85257871;v=rs509717;vdb=variation;vf=364091
Question: Hello everybody,

can somebody explain me what is the difference between replicates and conditions in RNA-seq? I found those terms in TopHat, Cufflinks and CummeRbund, but I do not understand it very well :-(

Sorry for stupid question.

Paul.

Answer: These terms will be common in all of science and might best be illustrated with an example. Say we have 20 mice and we evenly divide them into two groups of 10. In the first group we administer some sort of drug via injections while in the second group we administer a placebo. We then extract RNA from some tissue and analyse it (presumably via RNAseq, but the method and material don't matter). The RNA from the individual mice in each group would be replicates. The groups themselves represent different conditions (drug vs. vehicle/placebo). The terminology is the same regardless of material and also holds outside of biology.
Question: I have the following numbers, but I can't figure out where I found them. I can't find them again, and it's driving me crazy!

Thanks for your help!

Studies: 210832
Samples: 6322688
Runs: 7675139

Answer: The first can be found [here](https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=studies), the second [here](https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=samples), and the third is...somewhere else. Can't help you there, I couldn't find it.
Question: I am trying to figure out whether a compound is a potential drug precursor using biopython. And I need to search the atoms in the compound within a certain range of a given atom in a given residue. But I don't know how to use the NeighborSearch in the Bio.PDB and after reading this website(https://www.biostars.org/p/1816/) just made me more confused....

Could anyone procide me with a specific example about how to use the NeighborSearch? Thank you very much!

Answer: Normally you should write down what have you tried so others understand the specifics of the problem. That being said, here is an concrete example for NeighborSearch. Suppose you have a PDB `1prot.pdb` in the current directory and want to find all residues at distance 5 angstroms from alpha carbon in residue position `311` in chain `A`:

import Bio.PDB

parser = Bio.PDB.PDBParser(QUIET=True) # QUIET=True avoids comments on errors in the pdb.

structures = parser.get_structure('1prot', '1prot.pdb')
structure = structures[0] # 'structures' may contain several proteins in this case only one.
target_atom = structure['A'][311]['CA']

atoms = Bio.PDB.Selection.unfold_entities(structure, 'A')
ns = Bio.PDB.NeighborSearch(atoms)

close_atoms = ns.search(target_atom.coord, 5)

`close_atoms` is a list containing all atoms at distance 5 from the coordinates of `target_atom`. Note that this includes `target_atom` itself. To eliminate it from the list do: `close_atoms.remove(target_atom)`.
Question: Are all [GO terms][1] manually curated or are they computationally curated or mixture of both?

How about GO from [MSigDB][2] ?

[1]: http://www.geneontology.org/
[2]: http://software.broadinstitute.org/gsea/msigdb/

Answer: Yes it is manually created. From their own docs:

> Shared vocabularies are an important step towards unifying biological databases, but additional work is still necessary as knowledge changes, updates lag behind, and individual curators evaluate data differently. The GO aims to serve as a platform where curators can agree on stating how and why a specific term is used, and how to consistently apply it, for example, to establish relationships between gene products.
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: Here I query the UCSC mysql anonymous server for the coordinate of the region between the CDS and the transcription sites (5' UTR, but you can extend this position to get a longer 'promoter' ) ( only for strand "+", for the reverse strand use cdsEnd and txEnd...). It builds an cURL query for the USC DAS server. This url is then piped into *sh* to get the genomic sequences.

mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -D hg19 -N \
-e 'select concat("curl http://genome.ucsc.edu/cgi-bin/das/hg19/dna?segment=",chrom,":",txStart+1,",",cdsStart+1) from knownGene where strand="+" and txStart!= cdsStart limit 10' |\
sh > result.concatenated.xml

result:

cttgccgtcagccttttctttgacctcttctttctgttcatgtgtatttg
ctgtctcttagcccagacttcccgtgtcctttccaccgggcctttgagag
gtcacagggtcttgatgctgtggtcttcatctgcaggtgtctgacttcca
gcaactgctggcctgtgccagggtgcaagctgagcactggagtggagttt
tcctgtggagaggagccatgcctagagtgggatgggccattgttcatctt
ctggcccctgttgtctgcatgtaacttaataccacaaccaggcatagggg
aaagattggaggaaaga

gggtctccctctgttgtccaaggctggagtgtagtagtgctatcgcagct
gactgcagcctcaaccttccaggctgaagcgatcctcccacctcaacctc
ccacgtggctgagactacaggtgcttgccactatgcccaactaacatttg
gaattttcgtatacgtggattccagaggggtgacagcgaaacgtgagtaa
(...)
Question: Hi,

I have used both [DAVID][1] and [g:profiler][2] to do functional annotation and pathway analysis of a set of genes that I have. Now, you might call this redundant because both effectively does the same thing, but I am doing this to validate my results using two different tools.

However, I do not know the appropriate way in which to display my results. A casual google search is not throwing up many results which I can use to reliably pick a way in which to display my results. Kindly provide some papers or material I can use to display the results I have in a decent manner. Thank you.

[1]: https://david.ncifcrf.gov/
[2]: https://biit.cs.ut.ee/gprofiler/

Answer: Hi,

There are several ways to display that type of data. Two simple ways are:

1) Tables with GO terms and p-values (at the very least)

2) Bar plots of the top GO terms ordered by p-value

More "prettier" ways include:

3) Append GO results to a heatmap, just as the example shown in the `ComplexHeatmap` R package shown [here][1]:

4) Circos plots (you can get really creative here)

[1]: http://zuguang.de/supplementary/ComplexHeatmap-supplementary1-4/supplS1_TCGA_OncoPrint/supplS1_OncoPrint.html
Question: I want to know whether detected variants are hotspot mutations in TCGA or not.

Though cBioportal provides hotspot site for each gene, it doesn't provide download function.

Are there sites that provide hotspot mutation list or table?

Thank you

Answer: I don't know of a database. But you can ask whether your specific mutations are recurrent by querying cBio over the web, using Python.

import urllib2

# get information on breast cancer "brca" mutations in gene PIK3CA
cBioUrl = "http://www.cbioportal.org/webservice.do?cmd=getMutationData&genetic_profile_id=brca_tcga_mutations&gene_list=PIK3CA"

# read cBio to find PIK3CA mutation info in breast cancer, which you can parse to find recurrent mutations
cBioFile = urllib2.urlopen(cBioUrl)
recurrenceData = cBioFile.read() # parse through the recurrenceData string (column "amino_acid_change")

The `recurrenceData` variable contains a list of samples with PIK3CA mutations, and the "amino_acid_change" column will tell you sites at which TCGA samples have mutations.

To scale this up, you would just make a loop where you replace the "PIK3CA" in cBioUrl with your gene of interest, and replace "brca" with your cancer type of interest. You can also query 19 common cancer types by looping through the following list in Python.

pancan19_cancerIDs = ["ov", "ucec", "gbm", "brca", "luad",
"prad", "skcm", "stad", "paad", "cesc",
"blca", "hnsc", "kirp", "lgg", "laml",
"kirc", "thca", "lusc", "coadread"]
Question: Is anyone aware of a tool which allows one to apply the concept of "mutation hotspots" to 3D protein structures? It is common to define mutation hot spots at the DNA level where we define the hotspot as a region of DNA sequence with recurrent mutations. But what if we want to know if some mutations are "geographically" co-located in the 3D protein structure (e.g., in the same pocket) but not necessarily "sequentially" co-located. Even without pre-defined concepts of "pockets" if a protein structure is known then it should be possible to define some concept of physical distance between residues. Any suggestions?

Answer: Some additional options:

- [Mutation 3D][1]
- [HotSpot3D][2]
- [HotMAPS][3]

[1]: http://www.ncbi.nlm.nih.gov/pubmed/26841357
[2]: https://github.com/ding-lab/hotspot3d
[3]: https://github.com/KarchinLab/HotMAPS
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: cols=c(rep("red",2),rep("blue",3))
with(rbind(A,B), plot(log2(V2/V5), log2(V8/V11), xlab="log2(FC)", ylab="log2(FC)", xlim=c(-5,5), ylim=c(-5,5), col=cols))
abline(v=0, lty=2); abline (h=0, lty=2)

You could also just add the colors as a column, which would make life easier if using ggplot2.