cannin/biostars_qa huggingface dataset:
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: Hi,

So it depends on what you want to do, personally I prefer Python and find it very easy to build html interfaces to databases very quickly (its very easy to get something working quickly):

For example pylons has a very simple way of building templates for webpages and it uses SQLAlchemy (an object-relational mapper) to access the database side.

http://pylonshq.com/

or

http://www.djangoproject.com/

Both have relatively good AJAX support. I know you may want to use Perl ... but I thought I'd suggest that there are some very easy Python frameworks available for this.
Question: The topGO package appears to have a function `showSigOfNodes`, which plots the GO graph and highlights the enriched terms. I have done my GO enrichment analysis with GOseq and that packages doesn't provide a plot function like that. Is there any way the function from topGO can be made to accept the data from GOseq?

If not, what else should be used? I know REVIGO and it's nice but it shows graphs based on semantic similarity, not the GO hierarchy.

Is there anything else available in R/bioconductor?

Answer: If you just want to view the hierarchy of the GO with provided GOids, you can have a try the combination of function `makeGOGraph` from package AnnotationDbi and `plotGOTermGraph` from GOstats.

---

UPDATE: 20170510

```r
goIDs <- c("GO:0051130","GO:0019912","GO:0005783")
color <- c("lightblue","red","yellow")
pv <- c(0.0001,0.03,1e-10)

## get results
pp <- getAmigoTree(goIDs=goIDs,color=color,filename="example",picType="png",pvalues=pv)
```

The graph can be get from: http://pan.baidu.com/s/1c2lEkxy
Question: I'm sorry if this is a repeated question, but I continue to have doubts.

I created a blastdb like this:

makeblastdb -in input.fasta -dbtype nucl -title test_DB -parse_seqids -out test_DB

But I can't understand how to add the taxid in order to have the same result has if I used the all nt database.

Can someone clarify me? Thanks

Answer: your command should look like:

makeblastdb -in input.fasta -dbtype nucl -title test_DB -parse_seqids -taxid_map taxidmapfile -out test_DB

The taxidmap file is a text file consisting of two columns.
You can download the taxonomy id information here:

ftp://ftp.ncbi.nih.gov/pub/taxonomy/accession2taxid/nucl_gb.accession2taxid.gz

You need to unpack it and you can make a taxidmap file by doing (something like):

```
sed '1d' nucl_gb.accession2taxid | awk '{print $2" "$3}' > taxidmapfile
```
Question: I have some older VCF files that don't have the contig length set in the VCF header. This means that Picard and some other tools that are very strict with the VCF spec won't accept them.

The contig entries in the header should be

##contig=
##contig=

but they are

##contig=
##contig=

I know that I can manually fix this by doing the following steps.

- unzipping the file
- extracting the header
- lookup the contig lenghts in a fasta.fa.fai file
- adding the lenght to the contigs records in the header with vim
- re-header with bcftools
- bgzip and tabix the the re-headered vcf file

In my hands this works but if you make a slight VIM copy paste error you will have spend a lot of time reheadering and bgzipping a large VCF file for nothing.

Therefore I would like to have a more robust automatic solution where I just give the VCF file and the reference genome file and the header is automatically fixed and a new bgzipped VCF file is written out.

Is there a tool that can add the contig lenghts to the VCF header and write out a new bgzipped VCF file?

Answer: ## Use bcftools and a reference fasta index

You can use the following [bcftools command][1] to add contig info into your VCF / BCF files:

bcftools reheader --fai ref.fa.fai -Oz input.vcf.gz > output.vcf.gz

It will also remove unused contigs and generally fix issues.

[1]: https://samtools.github.io/bcftools/bcftools.html#reheader
Question: Hi,

I am trying to parse MEME formatted motif using Bio.motifs object.

>>> list(Motif.parse(open("test.meme"),"MEME"))
...
AttributeError: type object 'Motif' has no attribute 'parse'

Does anyone know if there is a replacement function for this?

Thanks.
J.

Answer: [Looking at the docs][1], your usage looks off

>>> from Bio import motifs
>>> with open("Motif/alignace.out") as handle:
... for m in motifs.parse(handle, "AlignAce"):
... print(m.consensus)

[1]: http://biopython.org/DIST/docs/api/Bio.motifs-module.html

So motifs, not Motif?
Question: Hi everyone,

Sometimes I am asked to mention the most exciting moment or part of my research works, I often say that I was excited first time I know how to work with command lines in Linux. However I think might be there is a specific purpose behind asking this question. What would you reply if you asked so?

Thank you very much

Answer: I believe that is a common job interview question, in general, it is used to check with the candidate what are his/her expectations and what consider important in his/her career. It also depends on the background, a biologist entering into bioinformatics will be excited about CS stuff and vice-versa. So, there is no correct answer, just be sincere.

On my particular case, I am constantly excited but I can highlight:

- Learning the command line
- Learning to code
- The first time I use an HPC
- The first paper in a high-impact journal
Question: Hello, everyone, I am trying to reproduce the lncRNA get from the paper "Portraying breast cancer with lncRNAs" by Olivier et. al. Here I came to the reproducing step:

![enter image description here][1]

I don't know whether my code is doing right, can someone give me your comment?

# load RAW files
RAW_files = list.files(path = "./RAWsets", pattern = "_RAW.tar")
# RNA.list = list of RAW file after processed by Frozen Robust Multiarray Analysis
RNA.list = list()
# iterate through the .CEL files
for (RAWstr in RAW_files) {
print(RAWstr)
# extract those raw .tar files
RAWdirectory = paste("./RAWsets", substr(RAWstr, 1, nchar(RAWstr) - 8), sep = "/")
tardir = paste("./RAWsets", RAWstr, sep = "/")
untar(tardir, exdir = RAWdirectory, verbose = F)
# extract those extracted .CEL files
cels = list.files(RAWdirectory, pattern = "CEL.gz", full.names = T)
sapply(cels, gunzip)
# read all cel files to data.raw variable
data.raw = ReadAffy(filenames = substr(cels, 1, nchar(cels)-3), verbose = F, cdfname = "hgu133plus2cdf")
# frozen robust multiarray analysis
data.matrix = frma(data.raw)
# add to RNA.list as described above
tmp = list(data.matrix)
names(tmp) = substr(RAWstr, 1, nchar(RAWstr) - 8)
RNA.list = append(RNA.list, tmp)
}

### Combine into 1 matrix of all data sets
RNAdata = as.matrix(exprs(RNA.list[[1]]))
for (i in 2:length(RNA.list)) {
RNAdata = CombineMatrix(RNAdata, as.matrix(exprs(RNA.list[[i]])))
}
# normalize between arrays on a data of fRMA before
RNAdata_scaled = normalizeBetweenArrays(RNAdata)
# because samples names are GSMxxxxxx.CEL => erase .CEL extension
colnames(RNAdata_scaled) = substr(colnames(RNAdata_scaled), 1, 9)

# load data from table S1, S2, S3, S7
load("tables_set.RData")
# mutual index between total data and tableS1 which have clinical information
tableS1 = tableS1[match(colnames(RNAdata_scaled), tableS1[, 1]), ]
# compute batch data using tumor vs normal tissue as covariate
batch = as.numeric(lapply(tableS1[, 3], function(str) if (str == "normal") return(1) else return(2)))
# applied Combat algorithm in sva library with default parameters to adjust data for batch effects
RNAdata_scaled = ComBat(dat = RNAdata_scaled, batch = batch)

Your help is really appreciated!

[1]: https://image.ibb.co/h3Mvmk/Screenshot_from_2017_08_28_11_23_14.png

Answer: Hey, it seems that what you are setting as batch is in reality the covariate value (which is defined as the condition or "variable of interest" in the linear model). This is the code for applying ComBat if you specify the batches and the conditions that you have in your dataset

sample <- Pheno$Sample.Name
batch <- Pheno$Batch
condition <- Pheno$Covariate.1
pdata <- data.frame(sample, batch, condition)
modmatrix <- model.matrix(~as.factor(condition), data=pdata)
combat_data_matrix <- ComBat(dat=yourmatrix, batch=batch, mod=modmatrix)
Question: I have been asked to recommend introductory books and resources to R and Bioconductor. My problem is just, I never read a book to learn R or Bioconductor, so I have no experience with this and cannot recommend one. I am interested in mainly introductory books, possibly targeting various groups of readers (computer scientists, molecular biologists, (bio-)statisticians), any recommendation appreciated.

For example, I used the following resources:

The R-manuals, especially the R intro

There are also a lot of contributed documents there on the R web site, but I didn't
use.

If a package from Bioconductor interests me, I read the package vignette.

I read the Bioconductor mailing list, that helps to see what other people use.

I have the "Venebles, Ripley. S Programming" book, that is hardly introductory.

Which books did you find helpful or completely useless to learn R/Bioconductor? For example: R Programming for Bioinformatics looks promising, anybody read it?

Or do you share my reluctance towards R-books and prefer online resources?

Answer: "Statistics Using R with Biological Examples" by Kim Seefeld helped me a lot

also "Applied Statistics for Bioinformatics using R" by Wim P. Krijnen
Question: When extracting reads aligned to a certain region, what is the difference between using samtools view and bamtools filter region?

Here is the code and the counts of the number of reads recovered. I don't understand why the number of reads recovered are so different between the two approaches. Any help??

```
samtools index IV75_realigned_reads.bam
samtools view -hb IV75_realigned_reads.bam 3R:15,000,000-15,100,000 > Control_region.bam
bamtools count -in Control_region.bam
# 28034
bamtools filter -region 3R:15,000,000-15,100,000 -in IV75_realigned_reads.bam -out test.bam
bamtools count -in test.bam
# 7523775
```

Answer: as far as I can read in the manual of bamtools, region are defined with two dots and not using a hyphen. You should try (and don't play with fire, remove the commas):

3R:15000000..15100000

> REGION string format: A proper REGION string can be formatted like any of the following
examples (where 'chr1' is the name of a reference (not its ID)and '' is any valid integer position
within that reference.):
> -region chr1
> only alignments on (entire) reference 'chr1'
> -region chr1:500
> only alignments overlapping the region starting at chr1:500 and continuing to the end of chr1
> -region chr1:500..1000
> only alignments overlapping the region starting at chr1:500 and continuing to chr1:1000
> -region chr1:500..chr3:750
> only alignments overlapping the region starting at chr1:500 and continuing to chr3:750. T
Question: I have been asked to recommend introductory books and resources to R and Bioconductor. My problem is just, I never read a book to learn R or Bioconductor, so I have no experience with this and cannot recommend one. I am interested in mainly introductory books, possibly targeting various groups of readers (computer scientists, molecular biologists, (bio-)statisticians), any recommendation appreciated.

For example, I used the following resources:

The R-manuals, especially the R intro

There are also a lot of contributed documents there on the R web site, but I didn't
use.

If a package from Bioconductor interests me, I read the package vignette.

I read the Bioconductor mailing list, that helps to see what other people use.

I have the "Venebles, Ripley. S Programming" book, that is hardly introductory.

Which books did you find helpful or completely useless to learn R/Bioconductor? For example: R Programming for Bioinformatics looks promising, anybody read it?

Or do you share my reluctance towards R-books and prefer online resources?

Answer: I have R in a Nutshell on my desk and I use it at least weekly if not daily. Its well indexed with a lot of great examples and helps me save time by not searching online. The BioC chapter is decent, but quite short. I also bought Data mashups in R which I found to be a great purchase and a fun way to learn about GIS capabilities and the Yahoo API in addition to R.
Question: Dear all,

I was wondering if there is any tool to calculate PSI from the BAM files. Basically what I want is-

I have gtf file with all the exons and I want to calculate for each exon (lets take exon filled with grey for example) how many reads mapped to junction of that exon (a and b in figure) and how many reads skip that exon (c in figure). So if we have, let's say, 10k exons in the gtf, the output should contain 10k rows with a,b and c as columns.

I used [MATS][1] and [bam2ssj][2] but did not get what I was expecting. I have gtf file with approx. 300k exons but MATS gives the output only for 3k exons. bam2ssj is calculating the values for introns rather than exons.

Any suggestions?

Thanks in advance.

[1]: http://rnaseq-mats.sourceforge.net/user_guide.htm
[2]: https://github.com/pervouchine/bam2ssj

Answer: I just tried to write something for this: https://github.com/lindenb/jvarkit/wiki/Biostar103303

I don't have any data to test though.
Question: Hi,

Here is the `flagstat` output of my bam file:

37750740 + 352032 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
1965916 + 0 duplicates
20118207 + 184214 mapped (53.29% : 52.33%)
37750740 + 352032 paired in sequencing
18875370 + 176016 read1
18875370 + 176016 read2
15175066 + 138352 properly paired (40.20% : 39.30%)
16901306 + 153876 with itself and mate mapped
3216901 + 30338 singletons (8.52% : 8.62%)
1538902 + 13814 with mate mapped to a different chr
561025 + 5046 with mate mapped to a different chr (mapQ>=5)

extracted **properly-paired reads** with `samtools view -bf 0x2 1.bam > 1pp.bam` and again cross-checked with `flagstat`:

15175098 + 138352 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
930902 + 0 duplicates
15175066 + 138352 mapped (100.00% : 100.00%)
15175098 + 138352 paired in sequencing
7587550 + 69176 read1
7587548 + 69176 read2
15175066 + 138352 properly paired (100.00% : 100.00%)
15175066 + 138352 with itself and mate mapped
0 + 0 singletons (0.00% : 0.00%)
2 + 0 with mate mapped to a different chr
2 + 0 with mate mapped to a different chr (mapQ>=5)

How can I extract **properly-paired QC-passed reads** instead of extracting only properly-paired?

Kindly guide.

Thanks

Answer: This should do it:

samtools view -b -f 2 -F 524 1.bam > 1pp.bam

You appear to have unmapped entries with assigned chromosomes. Presumably you used bwa, since it's the only tool I know of that produces things like that.
Question: As the title says, I'd like to know if there's a set convention for matching dna/rna/aa letter symbols (that is ACGT for DNA, ACGU for RNA, and a lot more letters for AAs) with some properly defined colors (for example, DNA letters might be A: green, C: blue, G: yellow/black and T: red)? The idea is that such colors would then be used in the areas for better recognition of the characters. Thanks!

Answer: As far as I know, there are no standards for colors. Each program uses slightly different colors.

You might consider using the same colors as JalView: http://www.jalview.org/help/html/colourSchemes/

Here are the available options for nucleotides in Geneious:

And here are the available options for amino acids in Geneious:
Question: I am looking for interactions between a protein (ACTB) and calcium and I want to use the paxtoolsr package. Is this possible? What code is necessary?

Answer: Here is code that would get you interactions that include actin or calcium from Pathway Commons:

library(paxtoolsr)

# READ AND FILTER PATHWAY COMMONS SIF INTERACTION TABLE ----
# With the file PathwayCommons12.reactome.hgnc.txt.gz downloaded from
# https://www.pathwaycommons.org/archives/PC2/v12/ and uncompressed
# NOTE: Either download directly or use downloadPc2()

# Read file
dat <- readSifnx("PathwayCommons12.reactome.hgnc.txt")

# Pathway Commons uses CHEBI identifiers for small molecules, chemicals, drugs, metabolites, etc.
# IDs:
# * Calcium: "CHEBI:29108" from https://www.ebi.ac.uk/chebi/searchId.do?chebiId=CHEBI:29108
# * Actins from https://www.genenames.org/data/genegroup/#!/group/929

ca2 <- "CHEBI:29108"
actins <- c("ACTA1", "ACTA2", "ACTB", "ACTBL2", "ACTC1", "ACTG1", "ACTG2")

# Subset the SIF interaction table to those entries with either calcium or actin
# Use idsBothParticipants=TRUE if it must be both source and target nodes must come from ids parameter
sif <- filterSif(dat$edges, ids=c(ca2, actins), idsBothParticipants=FALSE)
sif_ids_both <- filterSif(dat$edges, ids=c(ca2, actins), idsBothParticipants=TRUE)

write.table(sif, "calcium_actin_interactions.txt", sep="\t", quote=FALSE, row.names=FALSE)

# > head(sif)
# # A tibble: 6 × 7
# PARTICIPANT_A INTERACTION_TYPE PARTICIPANT_B INTERACTION_DATA_SO… INTERACTION_PUB… PATHWAY_NAMES MEDIATOR_IDS
#
# 1 ACTA2 in-complex-with CALD1 Reactome NA Smooth Muscl… http://path…
# 2 ACTA2 in-complex-with ITGA1 Reactome NA Smooth Muscl… http://path…
# 3 ACTA2 in-complex-with ITGB5 Reactome NA Smooth Muscl… http://path…
# 4 ACTA2 in-complex-with LMOD1 Reactome NA Smooth Muscl… http://path…
# 5 ACTA2 in-complex-with MYH11 Reactome NA Smooth Muscl… http://path…
# 6 ACTA2 in-complex-with MYL10 Reactome NA Smooth Muscl… http://path…
Question: Hello,

I'm really new if the field and maybe my question is a little bit naive.

I'm trying to annotate my SNPs by using VEP. The reason to do this is to find the nonsense SNPs on each genome, the synonymous and non-Synonymous.
The organism that I'm working on is sheep.
The thing that is a little bit confusing to me is the tag system that Ensembl uses.
For example, the "synonymous_variant" is for the synonymous SNPs, but I'm not so sure about the non-synonymous and the nonsense. I'm taking the "coding_sequence_variant" and "stop_gained", respectively. Am I right?
Also, I cannot identify the CNVs. Is there any particular tag for this?

A second issue that I faced is that for some gene IDs there is no information about the name of the gene (symbol tag). Is there any way to use somehow a list with these IDs and find the names of these genes?

Thank you very much in advance and I'm really sorry for the questions "bombing".

Answer: The tag system is based on the [Sequence Ontology][1] (SO) consequence terms. Non-synonymous is not a SO term. This should be referred to as missense_variants according to SO. Check the SO definitions and a diagram showing the location of variants on the [Calculated consequence variants][2] page. The nonsense is known as stop_gained. If you annotate CNVs (larger insertions or deletions for example), you will have the same SO consequence terms. These are some of the consequences for [copy_number_variation according to SO][3] As for your second issue, I'd guess there is no gene name for the sheep gene, but you will have the Ensembl stable ID, e.g. ENSOARG00000005819. If you sent some examples, it'd be easier to help.

[1]: http://www.sequenceontology.org/
[2]: http://www.ensembl.org/info/genome/variation/predicted_data.html#consequences
[3]: http://www.sequenceontology.org/browser/current_svn/term/SO:0001742
Question: Hi all,
I'm running about 50 proteins of a particular species (sub_proteins.fa) against all proteins for that specie (proteome.fa) using hmmer. I have done this using phmmer with the following command:

phmmer --tblout results.table sub_proteins.fa proteome.fa

However, I noticed that in the resulting file (results.table) there are multiple lines where the same protein is being compared to itself (I guess because the 50 query proteins are also found in the proteome file). See below for an example:

# target name accession query name accession E-value score bias E-value score bias exp reg clu ov env dom rep inc description of target
#------------------- ---------- -------------------- ---------- --------- ------ ----- --------- ------ ----- --- --- --- --- --- --- --- ---
YLL_767 - YLL_767 - 8.1e-177 580.6 1.4 1e-175 582.3 1.4 1.0 1 0 0 1 1 1 1

is there anyway to prevent the hmmer programs from doing this? The reason I am concerned about this is because after this step I am going to be combining the query and significant targets together into a model that will then be used to search against other species proteomes...will having redundant proteins affect my results?

Any help is greatly appreciated!

Answer: There is no way to prevent the program from detecting identical matches. In most instances that's the exact purpose of these programs!

As to whether that will affect your downstream results, it depends on what you mean exactly by `I am going to be combining the query and significant targets together into a model`. If you mean that you plan to build a hidden Markov model from your query + identified sequences, then it doesn't matter that some sequences in your alignment are going to be identical. Model building procedure automatically down-weighs identical (and even very similar) sequences, so the net effect will be as if you did not have duplicates.

If you want to convince yourself that this is a case, make an alignment of two identical sequences - or 10 identical sequences, for that matter. When you run `hmmbuild` on that alignment, it will print a summary saying how many sequences are in the alignment (`nseq`), followed by how many effective sequences (`eff_nseq`) are there. No matter how many identical sequences you put in your alignment, the `eff_nseq` will be 1 at most, but more likely smaller than 1. That means even though the alignment has X number of sequences, the HMM building program will act as if it has seen only one of them. Even for large alignments of non-identical sequences, say 2000 of them, `eff_nseq` will rarely go into double digits unless your sequences are truly diverse.
Question: In Vcf file, generate by samtools mpileup, there is an AD field (Allelic depth). This field is suppose to replace the previous DP4 field (Number of high-quality ref-forward, ref-reverse, alt-forward and alt-reverse bases).

But when I use -t DP4, AD in mpileup, I have a different number of values in the 2 fields:

GT:PL:DP:DP4:AD 0/1:27,0,107:5:4,0,1,0:4,1,0

And I don't understand why one data is missing in the AD field... Is it not supposed to be like the DP4, so:

> high-quality ref-forward, ref-reverse, alt-forward and alt-reverse bases

?

After that I am not sure: what there are 3 values in the AD field ? And what are the differences between AD and DP4?

Thanks for your help.

Rob

Answer: I think DP4 shows ref and all alt reads, but AD further separates the different alt alleles (a single record can have multiple alt alleles).

If you check the AD and DP4 descriptions in the VCF header, it should explain both fields.

I actually don't have AD in my samtools 1.3 VCF. They may have taken it out or it might depend on some other setting. Thus, I can't really confirm.