cannin/biostars_qa huggingface dataset:
Question: Hi everyone,

I have a list of snps with string chromosome names that I want to replace with a sequential numerical value (preparing a qqman input file). For example, turn this:

chrA
chrA
chrA
chrB
chrB
chrC
chrC
chrC
chrC
...
chrZ
chrZ

into this

1
1
1
2
2
3
3
3
3
..
26
26

At the moment I am doing this with a script, using something similar to this:

sed 's/chrA/1/g'
sed 's/chrB/2/g'
sed 's/chrC/3/g'
..
sed 's/chrZ/26/g'

However, my actual genome has a large number of contigs, so this command takes a lot of time to run. Is there a more efficient way to do this?

Answer: Here's a solution with `replace` subcommand ([usage](http://bioinf.shenwei.me/csvtk/usage/#replace))
of [csvtk](http://bioinf.shenwei.me/csvtk/download),
just download the `.tar.gz` file, decompress and you can run :)

First you have to prepare a mapping file, which is a **plain tab-separated text file**,
you can easily use a spreadsheet software to create and export.

$ more mapping.tsv
chrA 1
chrB 2
chrC 3

I guess the SNP data file should be a tab-delimited file too. Here's a dummy one:

$ more data.tsv
chrA A z
chrA A x
chrA A c
chrB B v
chrB B d
chrC C tx
chrC C t
chrC C x
chrC C z

Then use `csvtk` to edit the SNP data file:

$ ./csvtk -H -t replace -f 1 -p '(.+)' -r '{kv}' -k mapping.tsv data.tsv
[INFO] read key-value file: mapping.tsv
[INFO] 3 pairs of key-value loaded
1 A z
1 A x
1 A c
2 B v
2 B d
3 C tx
3 C t
3 C x
3 C z

The long-option version would be easier to understand:

./csvtk --no-header-row --tabs replace --fields 1 --pattern '(.+)' --replacement '{kv}' --kv-file mapping.tsv data.tsv

**PS: this is a general method not limited to this case. `sed` is good for single replacement, `csvtk replace -k` can handle multiple replacements well, which is written in [Go](https://golang.org/) with good performance.**

**PS2: [seqkit](http://bioinf.shenwei.me/seqkit/) has exactly the same function to handle FASTA/Q files.**
Question: Do you know any public scientific SQL server ?

for example, I would cite:

UCSC http://genome.ucsc.edu/FAQ/FAQdownloads#download29

ENSEMBL http://uswest.ensembl.org/info/data/mysql.html

GO http://www.geneontology.org/GO.database.shtml#mirrors

(I'll give a +1 to each correct answer)

Answer: 1000 Genomes: since June 16, 2011: http://www.1000genomes.org/public-ensembl-mysql-instance

mysql -h mysql-db.1000genomes.org -u anonymous -P 4272
Question: Hi All,

I did a transcriptome assembly of Illumina SE and PE reads using Trinity, but the N50 values for my assembly are very low. Here is the summary:

```
Total trinity 'genes': 748144
Total trinity transcripts: 916206
Percent GC: 43.83

Contig N10: 1765
Contig N20: 1064
Contig N30: 677
Contig N40: 477
Contig N50: 369

Median contig length: 253
Average contig: 364.39
Total assembled bases: 333854406
```

In addition to Illumina PE and SE data, I also have PacBio full-length transcripts data. Is there any way I can perform a hybrid transcriptome assembly? Need your ideas on it.

PBcR with Celera assembler seems to be one option that corrects PacBio reads and performs hybrid assembly, but 1) it seems to be significantly slow as mentioned by authors, 2) it is not clear from its documentation whether it will use Illumina reads beyond the reads correction step i.e. for final assembly and 3) no example provided by authors on how to use it for 'transcriptome' hybrid assembly. These issues with PBcR makes me think that it is a slow error correction tool which uses only the PacBio reads for final assembly. Are there any other options for de novo hybrid assembly of transcriptome?

Thanks
Bade

Answer: If your PacBio reads are already corrected, you can use them with Trinity (you don't mention if they are or not corrected). Quoting from this post from Trinity mailing list:

"You can incorporate corrected pacbio reads into Trinity using the --long_reads parameter."

rnaSPAdes is based on SPAdes, so it should be able to take both types of data, and MIRA as well, but both of them will probably choke if the transcriptome is too large or complex.

If your PacBio are not corrected, correct them and try one of the above suggestions.
Question: I have a large set of RNA-seq expression data in a DGEList object, and I want to plot the epxression data between two factors for specific genes temporally.

I started out by subsetting the data into a smaller matrix and then realised that was silly, and I should be able to plot it from the DGEList object that the data is stored in. Each timepoint has three replicates so I would also be looking to take a mean of those replicates before plotting. Would subsetting the data first still be the best option or am I missing a far quicker and easier option.

DGEList Count:

Gene Symbol Sample1 Sample2 Sample3 Sample4 Sample5 Sample6 Sample7 Sample8 ...
Gene1 54 55 53 78 79 74 81 82
Gene2 23 21 22 45 44 47 61 62
Gene3 74 75 73 81 82 80 83 88
Gene4 2 3 1 10 9 8 12 11
...

MetaData:

Sample Name ... Day
[1,] Sample1 D0
[2,] Sample2 D0
[3,] Sample3 D0
[4,] Sample4 D3
[5,] Sample5 D3
[6,] Sample6 D3
[7,] Sample7 D7
[8,] Sample8 D7
...

Using the examples above, what I am trying to do is draw an expression line plot for Gene 2 and Gene 3, including averaging the expression levels on each day - but as I have two factors the samples come from two factors and so would need to be separate,

Answer: Okay, I get the feeling that this does not have to be anything special for now (in terms of a 'polished' plot). So, you could try this:

df
time gene1 gene2 gene3
sample1 day1 1 2 3
sample2 day1 4 10 3
sample3 day2 1 2 3
sample4 day2 1 2 3
sample5 day3 1 2 3
sample6 day3 1 2 3

Summarise by mean:

df <- aggregate(df[,2:ncol(df)], df[1], mean)
df
time gene1 gene2 gene3
1 day1 2.5 6 3
2 day2 1.0 2 3
3 day3 1.0 2 3

Plot

plot(1, type="n", ylab="Expression", xlab="Day (1, 2, 3)", xlim=c(1,3), ylim=c(0,10))
lines(gene1 ~ time, data=df, lwd=2, col="royalblue")
lines(gene2 ~ time, data=df, lwd=2, col="red4")
lines(gene3 ~ time, data=df, lwd=2, col="forestgreen")
Question: I'm doing an analysis of sequence features in 18S rRNA sequences. I have downloaded data from SILVA database. Unfortunately, typical vertebrate has more than a single 18S sequence. I know that most of them are considered pseudogenes, but unless I check them one by one, there's no way to guess which one is a canonical (reference) version. I have roughly 1000 sequences to verify from 150 organisms, so I'd prefer not to do it by hand.

Where to find or how to identify canonical version of 18S rRNA?

Answer: The SSU gene of S. cerevisiae (Z75578) is usually used as canonical to mark start and end positions of V1-V9 regions.
http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0087624 there's table 3 in supplemental materials you need.
Question: For Tumor Mutation Burden (TMB) usually people calculates non_synonym_mutations / Mb.

In this Ilumina's tech report (https://www.illumina.com/content/dam/illumina-marketing/documents/products/whitepapers/trusight-tumor-170-tmb-analysis-white-paper-1170-2017-001.pdf) they calculate TMB for TruSight Tumor 170 with both, all mutations and non_synonym. What strikes me is that the all_mutations have a higher correlation with WES. Is this an artifact because the low number of genes and the fact that this is an in silico selection of the WES? other papers have good correlations with 300 and 500 genes but probably 170 is too low.

Results of In Silico Studies
WES data from 5336 TCGA samples were filtered and analyzed in
silico using TruSight Tumor 170. TMB estimated from the TruSight
Tumor 170 targeted regions showed a high correlation to TMB
estimated from WES, with R
2
correlation values of 0.91 for total
mutations (Figure 2A) and 0.90 for nonsynonymous mutations (Figure
2B)

Seems that if you want a proxy for neo-antigens formation for improving immunotherapies, then counting non_synonyms would give you a better proxy.

Doesn't matter if you use all_mut or synonym? Could someone point me to a good paper measuring these correlations? Which is the limit for the number of genes where correlation starts to drop?

Answer: Correlation is a rather poor evaluation metric this case, since it can be driven by hyper-mutated outliers that are easy to get right with small panels, i.e. the correlation depends on both the mutation rate and the panel size.

Including silent mutations will give you a more accurate total somatic mutation rate because you add a few more data points and reduce the sampling variance. This can - in theory - also provide a slightly more accurate estimate of the true mis-sense rate.

I don't think there is a paper, because it's straightforward to do the power calculation or the downsampling from TCGA data. Maybe the recent one from the Van Allen group at DFCI where they compare whole-exome and small and large panels is helpful.
Question: Hi,

I'm looking for an easy way to retrieve all the genes in a list that are associated with a certain GO term, preferably using R/Bioconductor packages.
I'm not interested in under/overrepresentation or enrichment.

For instance, say I have a list of 1000 genes and I want to create a sublist with only the genes known to be involved in 'heart development'.

Thanks!

Answer: Using biomaRt within R:

library(biomaRt)
ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl") #uses human ensembl annotations
#gets gene symbol, transcript_id and go_id for all genes annotated with GO:0007507
gene.data <- getBM(attributes=c('hgnc_symbol', 'ensembl_transcript_id', 'go_id'),
filters = 'go_id', values = 'GO:0007507', mart = ensembl)
Question: Hi folks! I need your help with a ggplot2 representation.

I have a data set that has the following structure:

log2FoldChange Sequence_biotype Knockdown
-1.40 LTR A
-1.11 DNA B
-3.46 Protein A
-1.25 Protein C
1.03 DNA B
... ... ...

I am plotting the foldChanges as boxplots, one for each `Sequence_biotype`. I have the knockdown variable faceted.
What I'm trying to do is to do a one sample Wilcoxon test for each boxplot, comparing the log2FoldChange to 0 (to see if there is a significant change). That can be achieved with the following code:

`wilcox.test(x = data, mu = 0)`

when the data is grouped by Sequence_biotype and Knockdown.

My question is, how could I introduce the results of the wilcox test to the plot as labels or how could I compute it directly in the plot (I have been trying with `stat_compare_means(method = "wilcox", paired = FALSE)` from ggpubr package but all the pvalues are piling up to the same spot)

Thank you before hand, any help is appreciated!

Best,

Jordi

Answer: In case anyone faces the same problem, I will paste the code I used to get the plot.

# Subset the dataset into different knockdowns
A = com_sig.df[com_sig.df$Knockdown=="A",]
B = com_sig.df[com_sig.df$Knockdown=="B",]
C = com_sig.df[com_sig.df$Knockdown=="C",]

# Compute the test for the first knockdown and the position in the plot
stat.test_A = compare_means(log2FoldChange ~ 1, paired = FALSE, data = A, method = "wilcox.test",
group.by = "Sequence_biotype", mu = 0) %>%
mutate(y.position = 4.5)

# Add a label to identify the knockdown
stat.test_A$Knockdown = "A"

# Same for second knockdown
stat.test_B = compare_means(log2FoldChange ~ 1, paired = FALSE, data = B, method = "wilcox.test",
group.by = "Sequence_biotype", mu = 0) %>%
mutate(y.position = 4.5)
stat.test_B$Knockdown = "B"

# Same for the third one
stat.test_C = compare_means(log2FoldChange ~ 1, paired = FALSE, data = C, method = "wilcox.test",
group.by = "Sequence_biotype", mu = 0) %>%
mutate(y.position = 4.5)
stat.test_C$Knockdown = "C"

# Combine all 3 datasets into 1
stat.test = do.call(rbind, list(stat.test_A, stat.test_B, stat.test_C))

# change column name of the variable used to compute the wilcox test
names(stat.test)[names(stat.test)==".y."] = "log2FoldChange"

com = ggboxplot(data = com_sig.df, x = "Sequence_biotype", y = "log2FoldChange", fill = "Sequence_biotype",
facet.by = "Knockdown") +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
theme_bw() +
theme(legend.position = "none",
axis.title.x = element_text(hjust = .975, vjust = -.6),
plot.title = element_text(hjust = .5),
axis.text.x = element_blank(),
strip.background = element_blank(),
strip.text.x = element_blank()) +
geom_text(data = stat.test, aes(y=y.position,label = p.signif))
com

[Plot Generated][1]

[1]:https://ibb.co/hXY4rbT
Question: Dear All,

I have a data matrix that has 17 samples and over 800 genes that belong to ten different gene families. I want to show these gene families in heatmap by marking them.

I performed heatmap but I do not know how to show gene families in the heatmap graph.

Anyone knows how to that?

Answer: Yes, as per Sean, in *ComplexHeatmap*, you can segregate your heatmap into 'blocks' of different genes using the 'split' parameter. You can end up with nice heatmaps like this: https://www.biostars.org/p/273155/#273172

**Edit April 16, 2019: skip to the working example, here: https://www.biostars.org/p/286187/#286507**
Question: My expression data frame pretty small subset of it is as such

gene value_1 value_2
XLOC_000060 3.662330 0.3350140
XLOC_000074 2.568130 0.0426299

I been struggling with ggplot how to trasform the data so that I can plot it in ggplot i basically how do I define the factor etc.

saw few examples but I guess im still not able to replicate what they did with the sample data

I end up getting this error `Error in eval(expr, envir, enclos)` its a pretty fundamental error I suppose which I fail to figure out..

Any help would be highly appreciated

Answer: You have to *tidy* your data to plot it effectively with `ggplot2`. There are many ways to do this, below I show two (assuming you want to plot a heatmap). Note the difference in data rearrangement between `d` and `dd`. Also, highly recommend reading the [tidyr vignette][1] and the *Tidy data* section in the book *"R for data science"* by [Hadley Wickham][2] (creator of ggplot2, tidyr and reshape2). The online version of this book is available [here][3].

# create sample (untidy) dataset.
d <- read.table(header = TRUE, text =
"gene value_1 value_2
XLOC_000060 3.662330 0.3350140
XLOC_000074 2.568130 0.0426299")

d
gene value_1 value_2
1 XLOC_000060 3.66233 0.3350140
2 XLOC_000074 2.56813 0.0426299

# tidy it using tidyr package.
library(tidyr)
dd <- d %>% gather(sample, value, -gene)
dd
gene sample value
1 XLOC_000060 value_1 3.6623300
2 XLOC_000074 value_1 2.5681300
3 XLOC_000060 value_2 0.3350140
4 XLOC_000074 value_2 0.0426299

# plot heatmap.
ggplot(dd, aes(x = sample, y = gene, fill = value)) + geom_tile()

# tidy it using reshape2 package.
library(reshape2)
dd <- melt(d, variable.name = "sample")
ggplot(dd, aes(x = sample, y = gene, fill = value)) + geom_tile()

[1]: https://cran.r-project.org/web/packages/tidyr/vignettes/tidy-data.html
[2]: http://hadley.nz
[3]: http://r4ds.had.co.nz
Question: Hi. I'm using pheatmap to create heatmaps. I've made a heatmap and used the clustering option (method = correlation) and the resulting heatmap has a certain row order based on the clustering. I want to make another heatmap using pheatmap that then has the same row order as the first one. Is there a way to specifically set the order of a second pheatmap based on the first?

I've been trying using the tree_row$order argument on the pheatmap object that I create. I think I can use this to set the order of the other, but then I don't know how to get pheatmap to display/write to file the changed pheatmap object.

Here's an example of some code:

> # create a pheatmap
>hm_1<-pheatmap(mat1, cluster_rows=TRUE) hm_2<-pheatmap(mat1, cluster_rows=FALSE)
>
> # get the row order of another pheatmap and set the first to the second
>hm_2$tree_row$order <-hm_1$tree_row$order

This works, but then I have this hm_2 object. How do I get it to display the resulting map again without calling the pheatmap argument?

And is there a batter way to do this?

Answer: you can reorder your second data based on the order of the first heatmap (where clustering has been applied)
in my case this was useful to display non-transformed values while plotting the heatmap based on transformed values

> 1) create the hm pheatmap using transformed and normalized data 2)
> extract row and col orders from it as you suggested
>
> # capture order for rows from first hm row.order =
> ph$tree_row$order col.order = ph$tree_col$order
>
> 3) reorder the non-transformed data based on the order of the data in
> the heatmap extracted above values.to.show = raw[row.order, col.order]
>
> 4) plot new pheatmap using the same data but overlay numbers from the raw matrix instead of the transf+norm numbers
> ... display_numbers=values.to.show, ...

Hope this helps and can be recycled!
Question: Hi all,

I need to run a PCA analysis on a microarray dataset to estimate which are the most relevant features. How can I do it? Is it possible?

I tried with `scikit-learn` but I was unable to come with the relevant genes. I did it like this:

from sklearn.decomposition import PCA
import numpy as np

# X is the matrix transposed (n samples on the rows, m features on the columns) - it is represented as a numpy array
X = np.array([[-1, -2, 5, 1], [-3, -1, 1, 0], [-3, -2, 0, 2], [1, 1, 1, 3], [2, 1, 1, 4], [3, 2, 0, 5]])

# suppose we want the 2 most relevant features
nf = 2
pca = PCA(n_components=nf)
pca.fit(X)

X_proj = pca.transform(X)

With `X`

array([[-1, -2, 5, 1],
[-3, -1, 1, 0],
[-3, -2, 0, 2],
[ 1, 1, 1, 3],
[ 2, 1, 1, 4],
[ 3, 2, 0, 5]])

it returns `X_proj`

array([[-2.9999967 , 3.26498171],
[-3.53939268, -1.18864266],
[-2.77013188, -2.15637734],
[ 1.67612209, 0.03059917],
[ 2.87464655, 0.35674472],
[ 4.75875261, -0.30730559]])

How can I say which are the selected features? Is there another way to do it (also in R for example)?

Thanks

Answer: Hi,

be careful because I think you are making some confusion about the concept of "components" in PCA.

The PCA doesn't return you the most relevant "features" in the dataset; it rather returns two vectors (or more, depending on how many components you calculate), corresponding to the coordinates of each element on a plane that divides the datasets into two parts. For a quick explanation of how PCA works, you can read this document or this article on Nature Biotechnology.

If you want to determine how much each of your original variables contribute to the components, you can either look at the loadings of the pca (I can't install scikit right now, but it may be inside the pca object), or you can calculate the correlation of each variable with each of the components, like described here. This should give you an idea of the "relevant" features, the features that contribute the most to divide the data into two.

Another tool to do PCA with python is orange, which gives you both a graphical interface (if you want to use it), and a python library.
Question: I am analyzing microarray data generated using Illumina Human HT 12 chips, and there were multiple batches as the samples were analyzed. The data I have has been through the 'standard' genome studio normalization steps, but has not been adjusted for any batch effects.

In analyses testing an outcome of interested against the expression values it is common to 'adjust' (include as an independent variable) for the batch effects using a factor variable.

I have also seen elsewhere that analysts may adjust for the relative log expression (RLE) mean to account for technical bias. RLE means are more commonly used to assess the batch effects using boxplots - I can see from boxplots in my data the a couple of the batches have significantly higher RLE means, bot not all.

My question is which method most accurately accounts for the technical variability introduced by the batches?

My feeling is that using the RLE mean values is best because, not only is this a linear variable, but it is actually based on the data! The batches may not necessarily have affected the expression, but to include them as covariates anyway must introduce some noise to the model. Whereas including the RLE mean values as a covariate, which are based exclusively on the expression data itself, will only account for the observed technical variation. Is this rationale logical? Have I overlooked anything? Many thanks.

Answer: This article compared several methods of removing batch effects:
http://www.plosone.org/article/info%3Adoi%2F10.1371%2Fjournal.pone.0017238

By their metrics, ComBat ( http://www.bu.edu/jlab/wp-assets/ComBat/Abstract.html ) performed the best.

ComBat is available in the R/Bioconductor package SVA: http://www.bioconductor.org/packages/release/bioc/html/sva.html

The advantage of using these methods over simpler approaches is that (some of them) they can remove batch effects while "protecting" your model
of interest.
Question: I have to generate a "toy dataset" with one human chromosome and short reads mapping on these contigs.

Which one would you pick? I would go for chromosome 21 since: 1) it is short, so less data 2) it has no gender bias (depth) like X/Y chromosome.

Is this a good pick or would you advise to go for another one?

Are there pre-processed datasets anywhere with contigs and corresponding reads filtered for just the specific chromosome?

Answer: Depends. What sort of data? What's the goal? Chromosome 21 is pretty sparse, gene-wise, even taking its size into account. Chr19 is usually my go to, as it's still quite short while having a decent gene density. Really depends on what your goal is. Avoiding X/Y is probably a good idea, otherwise anything after chromosome 13 are all pretty small.
Question: Hi everyone

Can you please help me to extract SAM file from SRA?

I took the dataset from here http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM1208162

Downloaded SRA file.

Then did sra-dump on this (as it is said on the page of dataset that reads are already aligned).

But as a result I got very small file (~50 Mb). When I tried to convert it to sorted BAM file I got:

```
[bam_header_read] EOF marker is absent. The input is probably truncated.
[sam_header_line_parse] expected '@XY', got [@HD VN:1.3]
Hint: The header tags must be tab-separated.
[samopen] no @SQ lines in the header.
```

I tried to look at summary information of this SRA file here: http://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR951914

And didn't see any information about alignment

Tried to look at alignment information by command:

vdb-dump ./SRR951915.sra | grep "ALIGNMENT_COUNT"

got an error

```
vdb-dump.2.1.7 int: data bad version while constructing page map within virtual database module - VCursorCellData( col:PANEL at row #1074 ) failed
fastq-dump command led to an error: data bad version while constructing page map within virtual database module - failed SRR951914.sra
```

Any idea?

Answer: Hi, yesterday I got a good output from fastq-dump.

Your command line is wrong. To run fastq-dump correctly you should know whether your reads are single or paired. In your case, we have single reads, then the command line will be:

./fastq-dump --split-spot SRR951914.sra

For paired reads:

./fastq-dump --split-files *file*.sra