cannin/biostars_qa huggingface dataset:
Question: Hi I am using a GEO dataset with Chip Seq experiments and it has an uploaded bedgraph file and since I need to use bed files for my analysis, is there a way to convert a bedgrpah to bed , also I don't quite understand the difference between both.

Answer: [Bedgraph][1] is BED-like, but puts the score data into the fourth column and often contains headers.

To convert to BED, strip out the headers with `grep` and put in a dummy ID value with `awk`:

$ grep -Ev '^(browser|track|#)' myData.bedgraph \
| awk ' \
​ BEGIN { OFS = "\t"; } \
{ \
print $1, $2, $3, "id-"NR, $4; \
} \
' - \
> myData.bed

The five-column BED data now follow the convention of putting the ID into the fourth column, and the score value into the fifth column.

[1]: http://genome.ucsc.edu/goldenPath/help/bedgraph.html
Question: Hi all,

What is the difference between GRCh and hg? I have small RNA interaction data with genome like below and they have used hg19 for mapping the reads. I want to use GTF file in Ensemble which is based on GRCh. but I was wondering the hg and GRCh are the same? If not, where can I find the GTF file for the hg19?

chromosome start end interaction ID alignment score of target part strand
chr2 1.9E+08 1.9E+08 human08246 32 +
chr16 5.2E+07 5.2E+07 human04514 32 -
chr3 1.1E+08 1.1E+08 human11115 36 +
chrM	14360	14382	human08122	34	+

Answer: I believe the genomic content for the two is identical, except for the mitochondrial contig.

The contig names are also different. GRCh37 names them `chr1`, `chr2`,,`chr3`, etc, while hg19 just has `1`, `2`, `3`.

Thus you can use the same GTF file for both (excluding mitochondrial, of course) if you do a simple replace operation for the contig names.
Question: I would like to retrieve the abundance histogram of a graph that I am building from a set of sequences.

I have seen how I can iterate on the graph and get the distribution by myself. However, it seems, from the options, that the distribution is computed while building the graph. The `Histogram` class comes with all the bells and whistles, and I would like to use it. Am I right? Can I access it?

I have also perused your examples, but I cannot find what I am looking for:

- Example `kmer/kmer13.cpp` shows how to build the histogram for a Bank.
- Example `debruijn/debruijn26.cpp` shows how the get the abundance of a node.
- Example `storage/storage6.cpp` shows how to get the distribution from a file generated by DSK.

Any pointers?

Answer: Hi Matthias,

Here is a quick and dirty answer.

First extract the histogram from the h5 file:

gatb-core/bin/h5dump -y -d histogram/histogram mygraph.h5 | grep "^\ *[0-9]" | paste - - | sed s/\ *//g | sed s/\,// > histo_data

Then plot the distribution:

Rscript plot-kmer-histo.R histo_data

With `plot-kmer-histo.R` (from Claire Lemaitre)

```r
#!/usr/bin/Rscript

#usage : ./plot-kmer-histo.R dsk.histo
# or RScript plot-kmer-histo.R dsk.histo
args=commandArgs(TRUE)
if (length(args)<1) {
cat("Usage:\n")
cat("./plot-kmer-histo.R dsk.histo [output.png]\n")
cat("1 obligatory argument :\n - kmer histo file (as output by DSK)\n")
cat("optional arguments :\n - output file (default = input_file.png)\n")
quit()
}

hist.file=args[1]

tab=read.table(hist.file)[,1:2]

if (length(args)>1) {
png.file=args[2]
}else{
png.file=paste(hist.file,".png",sep="")
}
bitmap(png.file,"png256",width=7,height=6,res=300)
suppressWarnings (plot(tab$V1,tab$V2,type="h",log="y",xlab="kmer coverage",ylab="count",main=""))
d=dev.off()
```
Question: I'm using the bedtools genomecov command to find coverage stats for a human genome alignment:

bedtools genomecov -ibam foo.bam

The [bedtools documentation][1] says that you need to specify a .genome file so the tool knows how big the chromosomes are using the "-g filename.genome" option. I haven't done this, but the command still works.

How does it know what the chromosome sizes are? Is it pointing to a .genome file somewhere in the system that I'm not aware of or is this simply not necessary contrary to the documentation?

Thanks

[1]: http://bedtools.readthedocs.org/en/latest/content/overview.html

Answer: With a BAM file, this is unnecssary. Sorry for the confusing docs.
Question: Recently `getGEO` function of the package GEOquery suddenly started throwing the following error:

```r
gse=getGEO("GSE106977",GSEMatrix=T)
#https://ftp.ncbi.nlm.nih.gov/geo/series/GSE106nnn/GSE106977/matrix/
#OK
#Found 2 file(s)
#/geo/series/GSE106nnn/GSE106977/
#Error in download.file(sprintf("https://ftp.ncbi.nlm.nih.gov/geo/series/%s/%s/matrix/%s", :
#cannot open destfile 'C:\Users\Marti\AppData\Local\Temp\Rtmp8cEqMH//geo/series/GSE106nnn/GSE106977', reason 'No such file or directory'
```

After some debugging I found that there was an error in the `getAndParseGSEMatrices` hidden function
To overcome this issue I made a patch, that fix the issue.

First copy this code into Rstudio or wordpad

```r
getAndParseGSEMatrices=function (GEO, destdir, AnnotGPL, getGPL = TRUE)
{
GEO <- toupper(GEO)
stub = gsub("\\d{1,3}$", "nnn", GEO, perl = TRUE)
gdsurl <- "https://ftp.ncbi.nlm.nih.gov/geo/series/%s/%s/matrix/"
b = getDirListing(sprintf(gdsurl, stub, GEO))
b=b[-1]
message(sprintf("Found %d file(s)", length(b)))
ret <- list()
for (i in 1:length(b)) {
message(b[i])
destfile = file.path(destdir, b[i])
if (file.exists(destfile)) {
message(sprintf("Using locally cached version: %s",
destfile))
}
else {
download.file(sprintf("https://ftp.ncbi.nlm.nih.gov/geo/series/%s/%s/matrix/%s",
stub, GEO, b[i]), destfile = destfile, mode = "wb",
method = getOption("download.file.method.GEOquery"))
}
ret[[b[i]]] <- parseGSEMatrix(destfile, destdir = destdir,
AnnotGPL = AnnotGPL, getGPL = getGPL)$eset
}
return(ret)
}
environment(getAndParseGSEMatrices)<-asNamespace("GEOquery")
assignInNamespace("getAndParseGSEMatrices", getAndParseGSEMatrices, ns="GEOquery")
```

Save the file as `GEOpatch.R` in your working directory and then, when loading the `GEOquery` library, source the saved file:

```r
library(GEOquery)
source("GEOpatch.R")
```

Now it should get going...

```r
gse=getGEO("GSE106977",GSEMatrix=T)
#https://ftp.ncbi.nlm.nih.gov/geo/series/GSE106nnn/GSE106977/matrix/
#OK
#Found 1 file(s)
#GSE106977_series_matrix.txt.gz
#trying URL 'https://ftp.ncbi.nlm.nih.gov/geo/series/GSE106nnn/GSE106977/matrix/GSE106977_series_matrix.txt.gz'
#Content type 'application/x-gzip' length 32707196 bytes (31.2 MB)

sessionInfo()
R version 3.3.2 (2016-10-31)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 16299)

locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

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

other attached packages:
[1] GEOquery_2.40.0 Biobase_2.34.0 BiocGenerics_0.20.0
[4] BiocInstaller_1.24.0

loaded via a namespace (and not attached):
[1] httr_1.3.1 R6_2.2.2 tools_3.3.2 RCurl_1.95-4.10
[5] bitops_1.0-6 XML_3.98-1.9
```

Answer: Upgrading *R* and *GEOquery* will fix the issue noted here.
Question: I am looking for input on ways to put numbers on how much different factors/covariates contribute to gene expression patterns. For example, let's say we have RNA-seq (or microarray) data from male and female mice that have received two different treatments (or no treatment) and the samples have been prepped in two different batches. Now we would like a bar plot (or similar) of how much the different covariates (1) sex, (2) treatment and (3) batch effects contribute, respectively, to the gene expression.

Some ideas:

Use PCA on the expression matrix and calculate the correlation between PC1 scores and the covariate vectors.

Use some machine learning algorithm to try to discriminate between male/female (or treatments, or batches), and evaluate by cross-validation how well this is working. Use the prediction performance as a measure of how strongly the covariate is reflected in the expression profile.

Is there a way to find out this sort of summarized importance from the linear models fitted as part of limma/edgeR/DESeq workflows?

I'd be grateful for any thoughts.

Answer: As a first step, it might be useful to remove genes with count 0 in more than 90% of your samples. Then obtain normalized and log-transformed expression estimates (logFPKM; maybe edgeR and voom() can be of use?) for all your samples and genes. I assume you already have obtained them somehow. Then prepare your data in such a way that you have a data frame where a row represents 1 measurement and the columns geneId, sex, treatment, batch, logFpkm (maybe aggregate() or melt() functions can be useful to reshape your data). Then perform anova and look at the sums of squares and the p-values of the F-test. Most likely, geneId has the biggest influence on your logFpkms and hopefully treatment is a good second.

fit <- lm(logFpkm ~ geneId + sex + treatment + batch,data=yourData)
anova(fit)

When you have many samples and many genes you can also take a random sample of your genes (say 1000) to get an idea and then scale up to bigger sample sizes.

Edit: Brentp, you are right. Its maybe better to exclude the geneId from the linear model
Question: Hi to all

I have CSV file with 11 columns and rows in it.

I read it into R as dataframe.

file1 <- read.csv("example-copy.csv", head = TRUE, sep =",",stringsAsFactors=F)

I want to extract rows with Names IND1 and IND3 with all the columns data. I tried with following code but i am not getting expected data and getting NA's.

Please found link to my example data here [dropbox link to csv](https://www.dropbox.com/s/i133sytqq6m099w/example%20-%20Copy.csv?dl=0).

# csv file

INDIVIDUAL,pop,M1,M2,M3,M4,M5,M6,M7,M8,M9
IND1,Aline,TT,CC,GG,GG,GG,CC,CC,TT,AA
IND2,Bline,TT,CC,GG,GG,GG,CC,AA,TT,AA
IND3,Rline,TT,GG,AA,GG,GG,TT,CC,CC,AA
IND4,Inbredline,TT,CC,AA,GG,GG,CC,CC,TT,AA
IND5,Bline,TT,CC,GG,GG,GG,CC,CC,TT,AA

#expected results

INDIVIDUAL,pop,M1,M2,M3,M4,M5,M6,M7,M8,M9
IND1,Aline,TT,CC,GG,GG,GG,CC,CC,TT,AA
IND3,Rline,TT,GG,AA,GG,GG,TT,CC,CC,AA

and simple my code is here

file1 <- read.csv("reddy - Copy.csv", head = TRUE, sep =",",stringsAsFactors=F)
file2<-file1[c("IND1","IND3"),]
File2

I do not know where I went wrong in this code. Can any help me to solve this. any help in this regard will be highly appreciated

Thanks in advance

Answer: Although not a bioinformatics Q (purely R programming), I'll give a hint

There are various way to subset a data.frame in R. Look here and digest them all as they will be extremely useful for you in future

https://www.statmethods.net/management/subset.html

In particular, you may want to try

mySubset <- file1[ file1$INDIVIDUAL %in% c("IND1", "IND3"), ]
Question: For example, when you analyse Untreated condition vs Day1Treatment condition with only these samples in your design you will get some size factors.

If you have also set up a Day6Treatment condition in your design (but you don't use it yet), and want to to do the same comparison Untreated vs Day1Treatment, the size factors will change. (taking into account Day6)

I thought size factors was only relative to library size of its sample , so why do they change dependent on the design file even if you don't use a set of samples.

I was doing analyse Day1Treatment vs Untreated & Day6Treatment vs Untreated with two separate design files.
But now I am wondering if it's better to have one design with all the samples, and do the two comparison to get same sizefactors because at the end of the day you finish with different differential genes detected.

Answer: SizeFactors are calculated based on all the samples in your dds object. So if you have 1 large experiment its better to put everything together and then perform comparisons for example using contrasts, rather than making new dds objects with subsets of your original dataset.
Question: I performed denovo assembly with Trinity using reads from heart mouse RNaseq. Than I mapped the transcriptome back to the reference genome with Blat . I also used Kallisto to Count the transcript abbundance in each sample. But now I want to know what Trinity ID's are known already in the annotation and what their names, and what is not annotated. How can I do that?

Answer: **Now, for a correct answer**.

You are doing mice which is a really well annotated organism. Trinity will be imprecise. You would be way better using Cufflinks with the reference genome and annotation, if you are looking for novel isoforms or things like that it will find them for you. Are trying to achieve something in particular?

Cufflinks will also give you the correspondence between its gene names and the official ones.
Question: Hi all,

I have 27 soil WGS metagenome datasets and I am trying to assemble them into contigs that are at least 1000-2000 kb long. Each dataset on its own is 20-30 Gigabytes of paired-end read fastq files.

I first tried the [Ray Meta assembler][1] because it's supposed to run well in parallel. I was able to do that for most datasets but have gotten very short contigs (most are <500 kb). Then I found [this paper][2] that suggests it does better for low-complexity datasets.

I also took a look at [Concoct][3] and I think the strategy sounds like it makes sense, but the code on their github pages is woefully outdated and I'm not sure how much of it is still maintained. Also, it suggests combining all datasets into one and then trying to assemble it (a "coassembly") and using that for the downstream analysis, but that approach will be computationally challenging since my datasets are so large.

If anyone has any experience with assembling complex soil datasets, I'd love to do some brain storming, so please reach out!

Thanks!

[1]: https://genomebiology.biomedcentral.com/articles/10.1186/gb-2012-13-12-r122
[2]: http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0169662
[3]: http://www.nature.com/nmeth/journal/v11/n11/full/nmeth.3103.html

Answer: The best approach for contiguity is generally to coassemble if you have sufficient resources. Some assemblers (HipMer/MetaHipMer, Ray, and Omega/Disco) can distribute that to spread the memory use across multiple nodes... since you've tried Ray, you might try the others too and see if they give better results. For a single node, we've found Megahit gives the best results with the lowest resource consumption.

You can also try approaches such as binning (using e.g. Metabat) and then assembling just reads that map to each individual bin. Normalization, error-correction, and/or discarding low-depth reads can also improve assemblies. With Ray and Disco, both error-correction and merging paired reads prior to assembly increases continuity.

But in general it's not a solved problem, so you'll have to experiment a lot! Don't expect great continuity, though; complex metagenomes often yield an L50 (length) of 200bp or less.

Note that you may be able to bin the raw reads using a binning tool based on depth covariance if the 27 samples are different (different conditions, location, time, etc).
Question: Hi,

In DESeq(1) it was possible to compare between conditions like this:

res = nbinomTest( cds, "1", "5" )

Where 1 and 5 are were the conditions I wished to compare (and in total there were 6 conditions).

I have just started using DESeq2 but I cannot work out how to replicate this analysis.

I tried "res <- results(cds,"1","5")" but that doesn't work. Does using the follow specify what is being compared?

colDATA(cuffcdsnorepdeseq2)$condition <- factor(colData(cuffcdsnorepdeseq2)$condition,levels=c("1","5"))

The problem is when I use this and then "res <- results(cds)" I am unsure about exactly what is being compared.

How can I do a pariwise comparison in DESeq2?

Thanks,

Jon.

Answer: You'll want to use the `contrasts=` part of the `results()` functions. So something like:

res <- results(cds, contrasts=c("conditions","firstCondition","fifthCondition"))

This way, you can be confident in which comparisons you're making. I recall reading that the contrasts functionality is relatively new, so you may need to update if you haven't already done so recently.

Update: In fact it appears that the contrasts functionality is new in what is currently the devel branch (Bioconductor 2.13).
Question: I have downloaded a GWAS sumstats file and have a list of SNP locations.

The information for each SNP location is condensed in column 1 of the sumstats file in the following format:

1:100000012_G_T

- 1 = Chromosome
- 100000012 = SNP location
- T = Allele 1
- G = Allele 2

I'm trying to run an LD Score analysis on this sumstats file, but I first need merge these SNPs with SNPs in the hapmap3 SNP list that has SNP IDs in the rs number format (obvs incompatible with the info in the sumstats file).

I'm looking for the most efficient way to obtain rsIDs for all the SNPs I have in the sumstats file (munging the file itself is not an issue).

Is there an R/bash method/package/command that can do this (I've checked out bed2vcf in R but don't think I have enough information for this to work)?

Or, alternatively, is there a way to generate a VCF file from the information I have here? I've looked into `bcftools annotate` to generate my rs IDs. I have already downloaded hg19 dbSNP VCF file that contains all the rsIDs and associated info, but not sure how to generate the second VCF file (from my sumstats SNP ID info) to run this command.

Answer: Hello,

here is a way which can work. First create a python3 script called `sumstat2vcf.py` with this code:

import re
import sys

for line in sys.stdin.readlines():
data = re.split(':|_', line.strip())
print(data[0], data[1], ".", data[2], data[3], ".", ".", ".", sep="\t")

It expect a list of the location in the format you give above. So we can use `cut` on your sumstat file to get the first column. If the first line of the file decscribes the columns we have to use also `tail` to get rid of it. The result can be streamed directly to the python script to create a vcf file:

cut -f1 sumstat.txt|tail -n+2|python sumstat2vcf.py|sort -k1,1V -k2,2g > sumstat.vcf

For annotating I prefer [snpSift][1]. To snpSift we can also pipe directly:

cut -f1 sumstat.txt|tail -n+2|python sumstat2vcf.py|sort -k1,1V -k2,2g|java -jar SnpSift.jar annotate -id dbSNP.vcf.gz > sumstat_annotated.vcf

If you now just wants a list of rs IDs we can go on using `cut` and `grep`:

cut -f1 sumstat.txt|tail -n+2|python sumstat2vcf.py|sort -k1,1V -k2,2g|java -jar SnpSift.jar annotate -id dbSNP.vcf.gz|cut -f3|grep -e "^rs" > rsId_list.txt

Good luck :)

fin swimmer

**EDIT:**

Based on Alex' answer you can use `awk` instead of the python script. Just replace `python sumstat2vcf.py` with

awk '{n=split($0,a,/[:_]/); print a[1]"\t"a[2]"\t.\t"a[3]"\t"a[4]"\t.\t.\t."}'

**EDIT2:**

Tidied up the python code a little bit.
[1]: http://snpeff.sourceforge.net/SnpSift.html#annotate
Question: Hello,

I have a data set looking like this:

Gene Symbol	a1	a2	b1	b2
Prss12	475	347	342	367
2010008C14Rik	3	0	0	0
Gm10819	216	129	164	126
Gm2245	25	23	14	28

and I am loading into EdgeR with the following script:

#!/usr/bin/R
counts = read.table("gene_count.tsv",sep="\t",row.names=1,header=T)
group <- c("a1","a2","b1","b2")
gene_or_transcript <- "genes_"#for filenames so I know what I'm reading
countS = counts[,1:dim(counts)[2]]
library(edgeR)
#now we'll do the DGEList command
y = DGEList(counts=counts[,1:dim(counts)[2]] ,gene=row.names(counts))#counts=counts[,2:dim(counts)[2]] removes $Gene.Symbol, i.e. 1st column
o<-order(rowSums(y$counts),decreasing=TRUE)#remove duplicates
y<-y[o,]
d<-duplicated(y$genes)
y<-y[!d,]
y$samples$lib.size <- colSums(y$counts)
rownames(y$counts) <- t(y$gene)# exchange ID
y<-calcNormFactors(y)
data.frame(Sample=colnames(y),group)
design<- model.matrix(~0+group)
#
y<-estimateGLMCommonDisp(y,design, method = 'deviance')
y<-estimateGLMTrendedDisp(y,design)
y<-estimateGLMTagwiseDisp(y,design)

but this produces an error

"Error in return(NA, ntags) : multi-argument returns are not permitted
In addition: Warning messages:
1: In estimateGLMCommonDisp.default(y = y$counts, design = design, :
No residual df: setting dispersion to NA
2: In estimateGLMTrendedDisp.default(y = y$counts, design = design, :
No residual df: cannot estimate dispersion"

why is this error produced? how can I fix it? I don't see anything in the EdgeR manual about it.

thanks, -DEC

Answer: You only have 4 samples, you can't then have 4 groups. I presume you meant `group <- c("a", "a", "b", "b")`.
Question: I have RNA-seq data from a control and experimental condition for an organism that lacks a genome and I want to measure differential expression. I have already tried doing de novo assembly of the transcriptome separately on each sample but the contigs don't match, so I don't know how to compare them. Does anyone have any suggestions for how to assess differential expression with or without assembly of the transcriptome?

Answer: Co-assemble the control and experimental results to produce a reference transcriptome. Then map the reads from each experiment onto these reference contigs and perform differential expression analysis based on these mapping results.
Question: I have a table like below, it is a bed file of genome coordinate, I would like to keep only rows with numbers.

Input:

1	141009669	141009952
9	141016322	141016973
GL000195.1	81719	82468
GL000195.1	142613	142923
GL000220.1	119445	119746
HG115_PATCH 101957832	101958132
HG1308_PATCH 130205069	130205369
HG1308_PATCH 130205406	130205773
HG748_PATCH 77577953	77578264
X	200983	202660
y 205180	205702

output:

1	141009669	141009952
9	141016322	141016973
X	200983	202660
y	205180	205702

Thanks in advance!

Answer: $ sed -n '/^[0-9,X,Y]/Ip' test.txt
or

$ sed -n '/^[^A-WZ]/Ip' test.txt

1	141009669	141009952
9	141016322	141016973
X	200983	202660
y	205180	205702

with [tsv-utils][1]:

$ tsv-filter --iregex '1:^[^A-WZ]' test.txt

1	141009669	141009952
9	141016322	141016973
X	200983	202660
y	205180	205702

[1]: https://github.com/eBay/tsv-utils
Question: Is there an authoritative database to find which genes (SNPs) are already known for which disease? With no limit on the study method. Thanks! (New to biology)

Answer: Have a look at [DisGeNET][1]. It is a very well done and well supported database that aggregates and curates information from other sources. We are happy to use it in Open PHACTS (which also offers another way to access it).

[1]: http://www.disgenet.org