cannin/biostars_qa huggingface dataset:
Question: Hi everyone

I'm making heatmap by using r studio. My data have 4 rows, 1 with name and 3 with values. I export excel sheet by selecting numeric.

> library(readxl)
> drug_class <- read_excel("C:/Users/hondalab/OneDrive/Desktop/16s NTR/Hiseq/drug_class.xlsx",
+ col_types = c("text", "numeric", "numeric",
+ "numeric"))

when i put *heatmap* command but i got this error.

> heatmap(drug_class)
Error in heatmap(drug_class) : 'x' must be a numeric matrix

*print* command showed this;

> print(drug_class)
# A tibble: 29 x 4
`drug class` `NTR 1 Rain` `NTR 2 Dry` `NTR 2 Rain`

1 aminocoumarin 0.00233 0.0239 0.00393
2 aminoglycoside 0.00559 0.136 0.0725
3 antibacterial free fatty acids 0 0 0.0000264
4 benzalkonium chloride 0 0.000615 0.000123
5 carbapenem 0.00255 0.0321 0.0184
6 cephalosporin 0.00444 0.0644 0.0485
7 cephamycin 0.00240 0.0264 0.00830
8 diaminopyrimidine 0.00244 0.0324 0.0412
9 fluoroquinolone 0.00637 0.0528 0.0317
10 fosfomycin 0 0 0.000412
# ... with 19 more rows

Can anybody help to solve this error. Thanks!

Answer: heatmap says it requires a numeric matrix, so you probably have to convert your tibble to a matrix.

drug_class <- as.data.frame(drug_class)
rownames(drug_class) <- drug_class[, 1]
drug_class[, 1] <- NULL
drug_class <- as.matrix(drug_class)

Tidyverse ansywer

library("tidyverse")

drug class <- drug_class %>%
column_to_rownames("drug class") %>%
as.matrix
Question: Hi,

I'm just an student so I'm very new in the bioinformatics terms. I've read the terms mapping and annotation in articles and webpages but i'm not sure about the difference between them and I cannot find any nice info in google. At this moment I understand:

Mapping; Alignment against the whole organism genome. The output is the position of the sequence in the genome.

Annotation: Alignment against a specific protein/nucleotide/etc database. The output is the protein name, or molecule.

Is this correct?

In that case. It is not the same? I mean, if we get the protein or the molecule we know the position in the genome isn't?

My doubt is because I read in some miRNA pipeline that they do mapping and annotation in different process.

Regards and Thanks in Advance

Answer: Mapping is usually the process of aligning short reads to a reference sequence, where the reference can be a genome or transcriptome.

On the other hand, In the Annotation, you go much deeper by incorporating biological meaning to the sequences. In wikipedia, it says "**Genome annotation** is the process of attaching biological information to [sequences][1]."

In Genome annotation, you may map your genomic reads to validate your reference sequences or your transcriptomic reads to find the exons and the regions which do not code. But, along with mapping, you can also annotate the genome by predicting the genes using gene predictors, repeat regions, predict non-coding genes, segmental duplications, Variant and SNP analysis, Identifying the pesudogenes etc. All these also come under annotation. So, note that it is different from mapping.

Whereas in the **Transcriptome annotation**, again, after **mapping** the reads, you would still be interested in knowing the gene expression analysis, isoforms, pathway analysis, gene family analysis, functional annotation (If the transcriptome is denovo, then assigning the function to the transcript). Hope this helps you.

~Prakki Rama.

[1]: http://en.wikipedia.org/wiki/DNA_sequence
Question: I'm trying to decorate some kegg pathway maps with expression data using the [KEGG Mapper from the kegg website][1].
In the first text field I enter the kegg pathway of interest (zma00940) and then I upload a text file with my test data (see below). I click 'exec' but I get an error: 0 does not exist. I guess there is some problem with the format of my data dont see what it could be. Any help would be appreciated.

zma:100381820	20
zma:100273579	5
zma:103627433	1

[1]: http://www.kegg.jp/kegg/tool/map_pathway3.html

Answer: Apparently the KEGG color mapper requires a header row:

#zma	count
zma:100381820	20
zma:100273579	5
zma:103627433	1

When I add that to the file, it seems to work.
Question: I'm working with some ChIPSeq data from the ENCODE consortium (i.e. `wgEncodeBroadHistoneK562H3k79me2StdPk.broadPeak.gz`) but there is no strand information. I was reading [this][1] and they say that they use MACS (which returns strand information) and in [their ChIPSeq guideline paper][2] I couldn't find anything specific about the strand reported in narrow/broad Peak files....

The reason I'm asking is because I'd like to use Homer for motif discovery but it asks for strand information. it works even without it, but to be very precise I think I'd need those info too.

Does anybody have some input/ideas on this?

Thanks :)

[1]: http://genome.ucsc.edu/ENCODE/protocols/dataStandards/ChIP_DNase_FAIRE_DNAme_v2_2010.pdf
[2]: http://genome.cshlp.org/content/22/9/1813.full

Answer: I think that strand is ignored, particularly for histones, as it is misleading. Histones don't bind to one strand of the DNA, the whole of the DNA is wrapped around them twice. Yes, the reads will be of a single strand because that's how sequencing works, but they represent both strands.
Question: Hi

I have thousands(1000's) of fasta files in one directory, I want to replace all the fasta file headers with the same keyword **>Gast_superba** ??

suggestions.

Answer: At the risk of sounding redundant after Wayne's answer, I think you should be able to use the following command in Unix:

sed 's/>.*/>Gast_superba/' *.fasta
Question: Hello everyone

I tried to install ggplots packages of R in Ubuntu server, but there were some errors which displayed as follows:

Warning messages:
1:In install.packages("ggplot2"):
Installation of package 'plyr' had non-zero exit status

2: In install.packages("ggplot2"):
Installation of package 'reshape2' had non-zero exit status

3: In install.packages("ggplot2"):
Installation of package 'scales' had non-zero exit status

4: In install.packages("ggplot2"):
Installation of package 'ggplot2' had non-zero exit status

Then I tried to use ggplot2 through command:

library('ggplot2')

But terminal told me that

Error in library (ggplot2): there is no package called 'ggplot2'

I have found many information about this problem on internet, but no available ways find and error was error.

If you have some suggestion, please tell me.

Thanks for your help.

Answer: You need to install Rcpp first. I'm a bit surprised this wasn't done automagically.
Question: Hello everyone,

I've searched a whole lot on the web and there doesn't seem to be a popular strategy of tackling this question:

> Given the PI_HAT, Z0, Z1 and Z2 values between 2 individuals, is it possible to narrow down on the possible relationship between the two, in cases where it's not obvious (PO/FS/UN)?

I'll elaborate.

The Zx values are fraction of relatedness measures, as in fraction sharing 0 copies, fraction sharing 1 copy and fraction sharing both copies between 2 individuals. Totally unrelated people would have `Z0` ~ 1 and duplicated samples/identical twins would have `Z2` ~ 1. PO, FS, UN and to an extent HS relationships have well defined values for all Zx as well as PI_HAT

However, once you get to 3rd degree relationships and beyond, (`PI_HAT`~0.125), it gets murky. At that point, even the detailed breakdown with the Zx values don't help us figure out what the exact nature of the relationship could be. It could be first cousins, or it could be grandparent-grandchild. Both have `PI_HAT`~0.125.

In this case, I wonder if we have a comprehensive listing of all possible relationships for each `PI_HAT` fraction (0.125, 0.0625, etc) that I could then combine with clinical information (such as age) and then take a call on the actual relationship?

The closest I've gotten is [this table in Wikipedia][1] (and other sources, but Wikipedia is quite good, ironically):

Any help would be appreciated, thank you!

--
Ram

[1]: https://en.wikipedia.org/wiki/Coefficient_of_relationship#Human_relationships

Answer: The `PI_HAT`/`Z0`/`Z1`/`Z2` values produced by `plink --genome` are generally not accurate enough to justify that kind of relationship inference. I recommend looking at [KING](http://people.virginia.edu/~wc9c/KING/) or [PC-Relate](https://www.bioconductor.org/packages/release/bioc/html/GENESIS.html) for this purpose instead.
Question: Dears,

What is the difference between **lfc=log2(1.5)** and **lcf=1.5** ?

topTable <- topTable(fit, coef=1, number=Inf, adjust.method="BH", **lfc=1.5**)

Which is the correct one to use and why?

Another question is what cut-off to use???

I see many articles using |log2fold change (FC)|≥1.5 and others using |log2fold change (FC)|≥1.0? with a Q-value < 0.05.

If you can tell me a reference about it, I will be grateful

Regards,

Answer: > What is the difference between lfc=log2(1.5) and lcf=1.5 ?

If you are asking about which one to use in topTable, the correct is `topTable( ..., lfc = 1.5 )`, as the fold-changes reported (and used for thresholding) are already log2(FC). However, topTable should not be used, and topTreat is to be preferred - read the topTreat manual page with `?topTreat`.

As for "correct" threshold, my view is that there isn't one. It is of course necessary to control for false-positive discovery, hence the q-value of 0.05 (or lower). But for fold-changes, I guess it depends on if one is interested in all changes, or only big changes in expression.

One of the benefits of using log2(FC) is that the values are intuitive: a log2(FC) of one means one treatment has double the expression of the other; a log2(FC) of 2 means one treatment is 4 times more expressed than the other, and so on. However, using log2(FC) = 1.5 means one treatment is approximately 2.828428 times more expressed than the other - as for arbitrary thresholds go, you can't go much more arbitrary than this.
Question: Hello! everyone:

I'm new to linux, here I got a problem:

I have a file file1 like:

3
6
7
9
12

and file2 which is tab-delimited:

chr1 3052600 3052800 1 E3
chr1 3052800 3053000 2 E3
chr1 3059400 3059600 3 E3
chr1 3059600 3059800 4 E3
chr1 3059800 3060000 5 E3
chr1 3062600 3062800 6 E3
chr1 3101000 3101200 7 E3
chr1 3105000 3105200 8 E3
chr1 3105200 3105400 9 E3
chr1 3116800 3117000 10 E2
chr1 3117000 3117200 11 E2
chr1 3164800 3165000 12 E2

and I want to extract the lines in file2 which its 4-th column equal the number in file1 like below:

chr1 3059400 3059600 3 E3
chr1 3062600 3062800 6 E3
chr1 3101000 3101200 7 E3
chr1 3105200 3105400 9 E3
chr1 3164800 3165000 12 E2

I have spent several hours including wrote a very slow python script, and I searched for the oneline solution, but I got nothing!

awk -v FS="\t" 'NR==FNR{rows[$1]++;next}(substr($NF,1,length($NF)-1) in rows)' fiel1 file2

Thanks a lot for some suggestions!

Answer: I assume your file2 is TAB-delimited.

sed -i 's/^/\t/' file1.txt
sed -i 's/$/\t/' file1.txt
grep -fF file1.txt file2.txt > combined.txt
Question: I am using datasets available here (GFF and Gene sequences):

- https://plants.ensembl.org/Glycine_max/Info/Index

Could someone please let me know how I could cite the source of this data?

Answer: The genebuild was done by JGI. Further information and the reference are here: https://plants.ensembl.org/Glycine_max/Info/Annotation/#assembly
Question: Hello,

Hope someone can help me on this. Basically I'm mapping probeIDs from the chip hgu133a.db to entrezIDs and I found that out of 22283 probesetIDs, 1266 maps to multiple entrezIDs. within it there are microRNAs, histones, pseudogenes and other weird stuff :D. So question is...how do I make the call, if it is possible, to which is the most reliable mapping for each of the multimapping probeset?

Hope this question is pertinent to this community.

Thank you
Bruno
PhD student in Genomics
Roslin Institute, Scotland

Answer: A common approach is to ignore (remove) probesets that map to multiple genes. Another approach is to simply treat the probesets as "measurements", do your differential expression analysis, and if one of these multi-mapping probesets shows up as "important", dig into the annotation for that probeset specifically by hand.
Question: For a given cancer type in the [NIH Cancer Genome Atlas][1], I visit [the data portal][1] and download UNC RNASeqV2, level 3 expression data. Specifically, I grab files that end with the extension `*.rsem.genes.normalized_results`

Each file contains one line per gene, with the gene name and (I assume) its normalized FPKM expression value. I am assuming these data are normalized FPKM based on the filename and the [UNC RNASeqV2 protocol description][2] hosted on TCGA.

My questions are:

1. Are these expression data really measured in FPKM?
2. If they are, how should I convert from FPKM to TPM, for all the expression values for a given gene?

[1]: https://tcga-data.nci.nih.gov/tcga/
[2]: https://tcga-data.nci.nih.gov/tcgafiles/ftp_auth/distro_ftpusers/anonymous/tumor/acc/cgcc/unc.edu/illuminahiseq_rnaseqv2/rnaseqv2/unc.edu_ACC.IlluminaHiSeq_RNASeqV2.Level_3.1.0.0/DESCRIPTION.txt

Answer: At the end of this blog post, a simple formula is provided to compute TPM from FPKM:

TPMi=( FPKMi / sum(FPKMj ) * 10^6

edit: well, from the protocol you linked, and also from this wiki, the UNC V2 RNA-Seq Workflow uses MapSplice+RSEM, so I guess measures are already given as TPM - check here and here.
Question: I try integrate PAML in a pipeline in order to analyse a large number of alignment. But, I met a problem with the python interface to codeML. Indeed, while codeML works perfectly when I call it directly, I don't manage to make it work through the codeml.py interface.
I launch Python3.5 in the directory "/home/PAML" and all files are in this directory.

I launch codeML through those command lines in python3.5:

>>> from Bio.Phylo.PAML import codeml
>>> cml = codeml.Codeml()
>>> cml.alignment = "lysozymeLarge.nuc"
>>> cml.tree = "lysozymeLarge.trees"
>>> cml.out_file = "results.out"
>>> cml.working_dir = "/home/PAML"
>>> cml.set_options(noisy = 4)
>>> cml.set_options(verbose = 5)
>>> cml.set_options(NSsites=[0])
>>> cml.set_options(runmode = 0)
>>> cml.set_options(seqtype = 1)
>>> cml.set_options(omega = 0.4)
>>> cml.set_options(kappa = 2)
>>> cml.set_options(fix_alpha = 1)
>>> cml.set_options(CodonFreq = 2)
>>>
>>> cml.run()

And, I get this error:

Traceback (most recent call last):
File "", line 1, in
File "/usr/local/lib/python3.5/site-packages/Bio/Phylo/PAML/codeml.py", line 186, in run
Paml.run(self, ctl_file, verbose, command)
File "/usr/local/lib/python3.5/site-packages/Bio/Phylo/PAML/_paml.py", line 145, in run
stdout=subprocess.PIPE)
File "/usr/local/lib/python3.5/subprocess.py", line 560, in call
with Popen(*popenargs, **kwargs) as p:
File "/usr/local/lib/python3.5/subprocess.py", line 950, in __init__
restore_signals, start_new_session)
File "/usr/local/lib/python3.5/subprocess.py", line 1544, in _execute_child
raise child_exception_type(errno_num, err_msg)
FileNotFoundError: [Errno 2] No such file or directory: 'codeml'

It doesn't find codeML, but it is integrated in the PATH, and works perfectly anywhere through the terminal with only this command:

codeml codeml.ctl

Do you have any idea what could be the problem?

Answer: Does it work if you specify the full, absolute path to the codeml binary within Biopython? E.g.:

cml.run(command="/home/PAML/codeml")
Question: Im trying to make a 3D PCA but i m not sure if its really working or not i cannot label my sample which i intend to here is my code anyone can help or any suggestion how can i label my samples in my final plot.

d <- read.csv("~/analysis/gene_fpkm/PRIMARY_AUX_ATPASE.csv",header = TRUE)
names(d)
model <- prcomp(t(d[,-1]), scale=TRUE)
plotPCA <- function(x, nGroup) {
n <- ncol(x)
}
fit <- hclust(dist(x), method="complete")
groups <- cutree(fit, k=nGroup)
if(n == 3) { # 3d plot
plot3d(x, col=groups, type="s", size=1, axes=F)
axes3d(edges=c("x--", "y--", "z"), lwd=3, axes.len=2, labels = TRUE)
mtext3d(text)
grid3d("x")
grid3d("y")
grid3d("z")
} else { # 2d plot
maxes <- apply(abs(x), 2, max)
rangeX <- c(-maxes[1], maxes[1])
rangeY <- c(-maxes[2], maxes[2])
plot(x, col=groups, pch=19, xlab=colnames(x)[1], ylab=colnames(x)[2], xlim=rangeX, ylim=rangeY)
lines(c(0,0), rangeX*2)
lines(rangeY*2, c(0,0))
}
}
plotPCA(model$x[,1:3], 5)

this is the code im running i cannot label my sample in which comes in the final 3D pca .any suggestion or help would be highly appreciated

Answer: Here's a [tutorial I wrote][1] using the rgl package mentioned by Jean. It includes 3D labels, exporting and animating. Hope it helps!

[1]: http://www.2-bitbio.com/2017/04/animated-3d-pca-plots-in-r.html
Question: Is there any software available for simulating a bacterial population which contains similar strains?

Answer: See these papers:

Assessing the reliability of eBURST using simulated populations with known ancestry

http://bmcmicrobiol.biomedcentral.com/articles/10.1186/1471-2180-7-30

Hierarchical and Spatially Explicit Clustering of DNA Sequences with BAPS Software

http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3670731/

SimBac: simulation of whole bacterial genomes with homologous recombination

http://mgen.microbiologyresearch.org/content/journal/mgen/10.1099/mgen.0.000044

Modeling Bacterial Population Growth from Stochastic Single-Cell Dynamics

http://www.ncbi.nlm.nih.gov/pmc/articles/PMC4136116/

Also this one:

https://labs.wsu.edu/populationgenomics/resources/software-simulation-code/

There are many of them. The choice depends upon what exactly you need.
Question: I have a very basic question. I have created a local (mysql) database in my system. I would like to know what would be the best option for a front end. I had used html before and want to try out something new. I might also add some perl code in between for some computational work. Kindly help me in this regard..

Answer: [?]

- MySQL or PostgreSQL as Back end for the data
- (X)HTML for the Front end
- PHP to make the link between the Front End and the Back End
- Javascript to add interactivity to the Front End. I also use the Javascript library called jQuery that makes life easier to deal with Javascript
- SVG if I want for instance to draw interactive KEGG Pathways

[?]
Link to ["Professional PHP Programming"](http://www.amazon.com/Professional-PHP-Programming-Jesus-Castagnetto/dp/B0000B0T0W/ref=sr_1_1?ie=UTF8&s=books&qid=1277292463&sr=1-1)
Question: Hi all,

I made BLASTp analysis on a specific genome against large databases. Here is my command line for example:

> blastp -query SEQUENCE.fasta -db database -out blast_results -outfmt 6 -evalue 0.01 -max_target_seqs 5

As expected, my output files are very large and I would like to **select only the matches with a percentage of identical matches of 60 or higher**.

Here is an example of an output file with only the first few lines:

Prokka_00003	VFDB_prot_DB_VFG012837(gi:82524546)_(ipaH2.5)_invasion_plasmid_antigen,_fragment_[Mxi-Spa_TTSS_effectors_controlled_by_MxiE_(CVF465)]_[Shigella_dysenteriae_Sd197]	28.947	114	77	3	134	247	342	451	7.17e-05	44.7
Prokka_00007	VFDB_prot_DB_VFG042740(gi:15599489)_(pprA)_two-component_sensor_PprA_[type_IV_pili_(AI098)]_[Pseudomonas_aeruginosa_PAO1]	29.496	278	175	5	463	735	654	915	8.18e-20	94.7
Prokka_00007	VFDB_prot_DB_VFG038759(gb|YP_857329)_(fleS/flrB)_two-component_sensor_kinase_[ND_(AI144)]_[Aeromonas_hydrophila_subsp._hydrophila_ATCC_7966]	24.194	372	224	11	366	734	17	333	4.81e-19	89.4
Prokka_00007	VFDB_prot_DB_VFG038760(gb|YP_008044159)_(flrB)_flagellar-related_two-component_sensor_kinase_[Polar_flagella_(VF0473)]_[Aeromonas_hydrophila_ML09-119]	24.194	372	224	11	366	734	17	333	6.71e-19	89.0
Prokka_00007	VFDB_prot_DB_VFG038762(gi:330830289)_(flrB)_two_component_system_flagellar_sensor_histidine_kinase_FlrB_[Polar_flagella_(CVF786)]_[Aeromonas_veronii_B565]	26.141	241	163	5	497	734	105	333	7.11e-18	85.9
Prokka_00007	VFDB_prot_DB_VFG038761(gi:145298503)_(flrB)_two-component_system_flagellar_sensor_histidine_kinase_FlrB_[Polar_flagella_(CVF786)]_[Aeromonas_salmonicida_subsp._salmonicida_A449]	25.311	241	165	5	497	734	105	333	1.28e-17	85.1

Do you know how would I be able to do that?

Thank you so much for your much appreciated help!!

Have a nice day

Answer: According to the blast documentation, percent-identity is given in the third column of the standard output format "6". So you can simply sort or filter by the third column. In your example, there is no entry that satisfies your threshold.

awk -F "\t" '{if ($3 >=60) print $0}' blast-output.txt

Something like the above might work.
Question: I am trying to use TCGAbiolinks package to analyze illumina methylation array data, but it is saying to me that the server is down, this is happening since last Friday night. Anyone have some information about it??

Answer: I solved this problem by using `devtools::install_github('BioinformaticsFMRP/TCGAbiolinks')` to install this package instead of `biocLite("TCGAbiolinks")`. You can try this way.
Question: Hello,

I have been trying the last days to make my database for blast work but I really don't know how to figure out this problem.

Begining from the start of the problem, I want to make a blastp against a protein database that I have created.

Here is an example of my sequence identifiers and how they are formatted:

>TRA_G345655|locus=scaffold45|99209|101155|-|translate_table|standard

I try to make the blast database with the following command:

makeblastdb -in database.fa -title database -dbtype prot -out protdb -parse_seqids

It does not work and it gives me the following error:

No volumes were created because no sequences were found.

BLAST Database creation error: Defline lacks a proper ID around line 1

When I try it without '-parse_seqids', it works fine and it generated the three files with pin, psq and phr.

I put these 3 files with my database.fa fasta file in one folder called DB and then I try to run blastp with the following command:

blastp -query query.fa -db DB -out proteins_blastp_1e-30_table.txt -evalue 1e-30 -outfmt 6

or with:

blastp -query query.fa -db DB/database.fa -out proteins_blastp_1e-30_table.txt -evalue 1e-30 -outfmt 6

or with:

blastp -query query.fa -db database -out proteins_blastp_1e-30_table.txt -evalue 1e-30 -outfmt 6

But it gives me the following error:

BLAST Database error: No alias or index file found for protein database [database] in search path

How can I solve this problem or where is the origin of this problem?

Answer: The documentation on how to build proper sequence ids is [here][1]

[1]: http://www.ncbi.nlm.nih.gov/toolkit/doc/book/ch_demo/#ch_demo.T5

In short: your fasta defline doesn't contain a valid identifier and needs to be changed to be parsable. A valid identifier consists of two or three fields separated by `|`.

E.g. for a ref seq entry:

ref|NM_010450.1|

For your case most likely the sequence has a local identifier, then you should change your defline to start with:

lcl|TRA_G345655

The second error is caused by:

makeblastdb -in database.fa -title database -dbtype prot -out protdb -parse_seqids

**so you also need to change this**:

blastp -query query.fa -db protdb -out proteins_blastp_1e-30_table.txt -evalue 1e-30 -outfmt 6