cannin/biostars_qa huggingface dataset:
Question: Hi all,

I have expression level of trancripts and most of them contain a number after [.].
As example:
ENST00000352998.3
When I search in Ensembl database, instead of **ENST00000352998.3** it shows me **ENST00000352998.7**

Does this means that these two are two different transcripts ? or I can simply ignore numbers after [.] ?

Answer: Hi Jack,

The numbers after the '.' is versioning of the stable IDs. They represent minor changes to the underlying gene model, without having to re-assign a new stable ID. There is some more information in the Ensembl documentation pages:

http://www.ensembl.org/info/genome/stable_ids/index.html

Best wishes

Ben
Ensembl Helpdesk
Question: I want to retrieve the canonical sequence for some proteins from ensembl, based upon their UniProt IDs.

```r
library(biomaRt)

ensemblMart <- useMart(
biomart = "ensembl",
dataset = "hsapiens_gene_ensembl"
)

(result <- getBM(
attributes = c("uniprot_swissprot", "peptide"),
filters = "uniprot_swissprot",
values = "P31749",
mart = ensemblMart,
bmHeader = TRUE
))
```

This particular protein has [two isoforms][1], both of which are returned. I only care about the canonical form, so how can I distinguish the isotopes?

I suspect that I need to return another of the attributes from the mart

```r
listAttributes(ensemblMart)
```

but I don't know which one.

[1]: http://www.uniprot.org/uniprot/P31749#sequences

Answer: There is no canonical defined in BioMart.
Question: I've run FreeBayes on a few samples, but I'm not sure how to interpret the GL (genotype likelihood) part of the VCF format. As an example, here is the info from one of the samples:

GT:DP:AD:RO:QR:AO:QA:GL 0/1:28:13,15:13:521:15:611:-46.8976,0,-38.8067

I understand that there are 13 reads with the reference allele and 15 reads with the alternate allele, but I don't understand exactly what the p-values represented in the GL part are the significance values for, and why there are three values. The description I've found for GL is "Genotype Likelihood, log10-scaled likelihoods of the data given the called genotype for each possible genotype generated from the reference and alternate alleles given the sample ploidy.". So is it the p-values of 0/0, 1/1, and 0/1? What is the order? And why is one of them always log10 p-value of 0?

Answer: Ok, it seems like the order is 0/0, 0/1, 1/1, but the numbers are normalized by the most likely genotype, which will therefore always be 0.
Question: I am trying to convert plink files to VCF by using following commands

```
#!/bin/sh
# 1. have plink binary to specify reference allele
plink --noweb --bfile $plink_file --reference-allele $ref_Allele_file --make-bed --out $plink_file_modified

# 2. create plinkseq project
pseq $pseq_project new-project

#3. load plink file into plink/seq
pseq $pseq_project load-plink --file $plink_file_modified --id $plink_file_modified

#4. write out vcf file
pseq $pseq_project write-vcf | gzip > $plink_file_modified.vcf.gz
```

I have site-positions-file which looks like as follows

```
7 7672918
7 4458845
7 3013804
7 4518570
7 7216328
7 7523768
```

I need reference-allele-file like this

```
7 7672918 C
7 4458845 C
7 3013804 T
7 4518570 T
7 7216328 C
7 7523768 C
```

Kindly help me

Answer: There are several older threads on this, like

https://www.biostars.org/p/47561/

or

https://www.biostars.org/p/5862/

I guess the easiest way to circumvent your few steps and the use of pseq is to just download the PLINK 1.9 beta and use the recode function

plink --file your_ped_map_input --recode vcf
Question: I was a bit confused by the output of a bwa read-alignment. I have observations that return a good CIGAR value but poor MAPQ value. This alignment was obtained from converting a BAM file -> BED file that was generated with bwa mem.

chr1 142536052 142536087 336_1362_95 0 - 16 35M * 0 0 CACTCCAGACTCGGTGACAGAGTGAGATCCTGTCG 8?@?.:/@@@@@@@@@@@@@@@>?@@@@?@@>@=@ NM:i:1 MD:Z:8T26 AS:i:30 XS:i:29

I thought CIGAR values and MAPQ value should be correlated but this breaks my intuition. Any insights would be super helpful. Thanks!

Answer: The cigar string tells you how well the alignment matches the reference. That doesn't mean the nucleotide match, but that's not the point here.

The mapping quality tells you how likely this location is indeed the correct location for this read. A read matching perfectly to a repetitive element will have a low mapping quality because you are unsure about its exact location, it's a multi-mapping read.
Question: Dear all,

I'm trying to open a wiggle file from a ribosome profiling study with IGB (It's the first time I'm using IGB).

For example, the GEO accession is: [GSE53767][1].

I decompressed the ".wig.gz" file, but after trying to open it in IGB, it says:

Unable to load data set for this file.

After some searches, I found that a wiggle file should start with a declaration line like this:

variableStep chrom=chr1

This file didn't have it so I add this line, and then IGB was able to open it. But I wasn't able to load the data in IGB, it's like there was no data...

How can I properly open this file?

Thank you very much for your support.

[1]: http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE53767

Answer: I guess that the name of the chromosome you entered is not the one IGB has. look at the right hand side of the screen and you'll see the names of the available chromosomes, see if yours is similar to the the one you're viewing. I had the same problem.
Question: Hello.

I'm doing some microarray analysis, and I figured out that some probe ids have not any GENE SYMBOL. Am I suppose to delete those probes or let them exist ?

The way I apply annotation info for an Affymetrix HGU133plus2 array after RMA application, is the following:

probes=row.names(expressions)
Symbols = unlist(mget(probes, hgu133plus2SYMBOL, ifnotfound=NA))
Entrez_IDs = unlist(mget(probes, hgu133plus2ENTREZID, ifnotfound=NA))
expressions=cbind(probes,Symbols,Entrez_IDs,expressions)

Is something wrong with my code or such behavior is expected ? What do you guys do with these NA genes ?

Answer: It’s not unusual in older array. Due to changes in gene annotation and reference sequence, some probe set may be outdated (probes not mapping to any gene in the new build). To find a gene for your probe set, you need to blat all the probes of the probe set and look where they are mapping. ​
Question: I have paired end RNASeq Data of 8 samples (Illumina Hi - Seq). My question is whether I should trim off the adapters from these reads before to do alignment. They provided these adapters:

a) TruSeq Universal Adapter

5’ AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT 3'

b) TruSeq Adapter, Index Adapter

5’ GATCGGAAGAGCACACGTCTGAACTCCAGTCACXXXXXXATCTCGTATGCCGTCTTCTGCTTG 3'

If trimming is necessary, since these are paired end, does it mean I trim the Adapter-a from one of the pair (R1) and trim Adapter-b from another one of the pair? And which specific site 5' or 3' should be trimmed?

Additionally, these reads are 125bp reads, does it mean that `tophat` setting "`--mate-inner-dist`" should be specified as 125?

Thanks in advance!

Answer: The 3' ends of the reads should be trimmed.

Tophat mate-inner-dist is the distance bet the 2 reads of a pair, and depends on the average fragment length, you should obtain an idea of that value from the sequence provider, or whoever did the library prep.
Question: I am having trouble converting some MultiQC reports from HTML to PDF format while preserving the files' graphics. I am working on a Red Hat Linux cluster and trying to accomplish this with a Bash script. Currently, I have tried using several versions of wkhtmltopdf to do this, and, although it converts the files, the graphs and charts are not rendered. For example, the charts for where my Sequence Quality Histograms or Per Sequence Quality Scores should be will instead have a gray box with the text "loading..." inside. Has anyone found a tool that works for this?

Answer: The problem is that MultiQC plots are by default rendered as interactive JavaScript objects. These aren't traditional graphics and tend to not do well with HTML conversion tools such as the one you're trying. MultiQC can render static image plots though, using the option `--flat`. There's also a template called "simple" (use `-t simple`) which produces a report with no JavaScript at all.

The `--pdf` option renders a HTML report using the simple template and then passes this through the pandoc command (see the command [here](https://github.com/ewels/MultiQC/blob/44cf7ea9c8b153be1baebec76ebba0370aad8333/scripts/multiqc#L700-L709)). It's a little flaky as it requires you to install pandoc yourself, which can vary quite a bit between systems. It is something that could potentially be dropped in future versions of MultiQC.

Generally when people ask about PDF conversion I ask "Why?". I'm always curious - what's the use case here?

Phil _(author of MultiQC)_
Question: I have a set of data that looks like this:

```
NK.Chr1:75500000-95000000:28960-29007 NG-unitig0655 97.872 47 1 0 1 47 121009 120963 2.90e-14 80.6
NK.Chr1:75500000-95000000:28960-29007 NG-1DRT-unitig0549 97.872 47 1 0 1 47 623680 623726 2.90e-14 80.6
NK.Chr1:75500000-95000000:28960-29007 NG-1DRT-unitig0278 97.872 47 1 0 1 47 1224581 1224627 2.90e-14 80.6
NK.Chr1:75500000-95000000:28960-29007 NG-1DRT-Chr4 97.872 47 1 0 1 47 8416368 8416414 2.90e-14 80.6
NK.Chr1:75500000-95000000:28960-29007 NG-1DRT-Chr4 97.872 47 1 0 1 47 20041035 20041081 2.90e-14 80.6
NK.Chr1:75500000-95000000:28960-29007 NG-1DRT-Chr4 97.872 47 1 0 1 47 35175472 35175426 2.90e-14 80.6
NK.Chr1:75500000-95000000:28960-29007 NG-1DRT-Chr4 97.872 47 1 0 1 47 56460095 56460049 2.90e-14 80.6
```

I need to filter the lines in the range of 0-3900000, considering only the numbers before NG.

```grep 'NK.Chr1:75500000-95000000:[0-3900000]' NG.1DRT-blast.out > chr1-blast-NG.txt```

I tried this code, but it returned all the lines with ```NK.Chr1:75500000-95000000:```, not considering the range.

Anyone knows how to build a proper code for it?

Answer: I think the following should work:

awk -F '[:-]' '($1~/NK.Chr1/ && $2=75500000 && $3=95000000 && $4<=3900000 && $5<=3900000)'
Question: Hello !

I have a DESeq2 result dataframe-like structure in R with some NAs in the *p.adj* column. Strangely, those NAs are handled seemingly randomly by two extremely similar functions used to test the genes dow- or up- regulation. What is happening here ? A quick fix would be to change those NAs into "1" (never significant) but I want to understand :) Here is my code:

significant_D <- function(x){return(x$padj < 0.01 & (x$log2FoldChange) < -0.584)}

significant_O <- function(x){return(x$padj < 0.01 & (x$log2FoldChange) > 0.584)}

head(DESeq_results[which(is.na(DESeq_results$padj)),])
baseMean log2FoldChange lfcSE stat pvalue padj

WBGene00021406 20.718704 -0.21384520 0.5063939 -0.42229021 0.6728132 NA
WBGene00021407 3.961096 0.66807041 1.1159760 0.59864226 0.5494115 NA
WBGene00021405 1.939649 -1.16416923 1.6007395 -0.72726966 0.4670608 NA
WBGene00021409 1.719862 1.91952086 2.2842210 0.84033940 0.4007181 NA
WBGene00235257 7.055687 0.23653150 0.8488041 0.27866442 0.7805024 NA
WBGene00015246 15.699319 0.06366663 0.6514634 0.09772863 0.9221478 NA

sum(is.na(significant_D(DESeq_results)))
[1] 2558

sum(is.na(significant_O(DESeq_results)))
[1] 1452

sum(is.na(significant_D(DESeq_results)) & is.na(significant_O(DESeq_results)))
[1] 0

sum(is.na(DESeq_results$padj))
[1] 8582

Answer: The confusion is due to the following:

> TRUE & NA
NA
> FALSE & NA
FALSE
> NA & NA
NA

The output of `significant_D` and `significant_O` are boolean vectors with some `NA` values. It will necessarily be true that any `NA` or `TRUE` values output by one of these will be `FALSE` in the other (after all, the only time you can get a `TRUE` or `NA` is when the fold-change passes filtering, which means it will fail in the other function). Since `FALSE & NA` is `FALSE` the penultimate sum is 0. As an aside, it makes sense that `FALSE & NA` is `FALSE`, since `NA` can be considered "unknown".
Question: What does it mean gene sets are upregulated in phenotype for example cancer ?

After performing GSEA for phenotypes normal vs cancer, the report has been obtained and in it:

35 / 50 gene sets are upregulated in phenotype cancer. I know the meaning of a gene is up regulated or down regulated but what does it mean gene sets are up regulated? does it mean that all genes inside the gene sets were up regulated ? or what ?

Thanks

Answer: Hello

GSEA is a tool that determines whether a set of genes defined *a priori* shows concordant and statistically significant differences between two biological states. For that reason, upregulated gene sets mean that the genes contained in these gene sets are (the vast majority) upregulated, so the program determines that this gene set is upregulated.

Best
Question: In VCF files, what character encodings are used in practice?

What character encodings are *expected *for the format?

The examples I've seen seem to include single-byte characters. It's not clear whether it's intended to be 7-bit US-ASCII or 8-bit Latin1.

The [4.2 format specification][1] doesn't mention encodings for VCF files. (It does specify that Unicode characters are not supported in characters and strings in *BCF *files)

[1]: http://samtools.github.io/hts-specs/VCFv4.2.pdf

Answer: https://github.com/samtools/hts-specs/blob/VCFv4.3/VCFv4.3.tex

> The character encoding of VCF files is UTF-8. UTF-8 is a multi-byte character encoding that is a strict superset of 7-bit ASCII and has the property that none of the bytes in any multi-byte characters are 7-bit ASCII bytes. As a result, most software that processes VCF files does not have to be aware of the possible presence of multi-byte UTF-8 characters.
Question: I have multiple genes collected from diabetic patients, I need to find out the pathways related to these genes using KEGG, Reactome and Wikipathways.
I know how to do that but need to find out if the resulted pathways are related to Diabetes or not.
Please help me to do that, is there a tool or R package help me to find out diabetes pathways ?

Answer: The [rWikiPathways package in bioconductor][1] lets you search the database in a few different ways. Like the comments above, exactly how *you* define a "diabetes pathway" is up to you. For example, you can search for any pathway that mentions "diabetes" anywhere in the title or description and then filter for the human set:

library(rWikiPathways)
db.pathways <- findPathwaysByText('diabetes')
human.db.pathways <- db.pathways %>%
dplyr::filter(species == "Homo sapiens")
human.db.pathways$name

There were 4 that popped up when I ran this snippet:

> [1] "Type II diabetes mellitus"

> [2] "Insulin signalling in human adipocytes (diabetic condition)"

> [3] "Polyol Pathway"

>[4] "16p11.2 distal deletion syndrome"

You can retrieve the WPID identifiers and other info from the same returned data frame.

Alternatively, you can seach by your favorite diabetes gene, e.g., INSR:

findPathwaysByXref('INSR','H') # H denotes HGNC identifer for query xref

This returns 27 human pathways! The [bioconductor package][2] includes additional vignettes for other common use cases.

[1]: https://bioconductor.org/packages/release/bioc/html/rWikiPathways.html
[2]: https://bioconductor.org/packages/release/bioc/html/rWikiPathways.html
Question: Hello all,

**Introduction**

I am a recent graduate working for a Public Health Laboratory. I'm relatively new to bioinformatics, and most of what I know is based around NGS analysis. My lab director loves to challenge me. He wants to know different ways NGS (we have a MiSeq) could be implemented in our lab as an outbreak investigation tool.

I am aiming to do a study of Carbapenem-Resistant Enterobacteriaceae (CRE). The main goal is to be able to receive a CRE sample and use NGS to detect the genes (beta-lactamases) that are responsible. The idea is to be able to run quick analysis, while also compiling genetic information that could be used to connect the dots in an outbreak investigation (Phylogeny).

**What I already Know**

- The genes that I am looking for can be found within the bacterial chromosome, or within its plasmids
- For each gene, primers need to be designed for them.

**The Actual Questions**

First:
If I wish to take the whole-genome-sequencing approach, how would I be able to tell which parts of my output (FastQ) are plasmids vs. which parts are chromosome?

Second:
If I didn't want to do whole-genome, would it be possible to only sequence the genes that I'm looking for (if they are there)? And if so, how would I do it?

**Open for Discussion**

If anyone has any suggestions, solutions, or wishes to point me in a direction where I can learn more, please let me know. It would be a huge help, and is deeply appreciated.

Edit: ***Solution Found***

Thanks to those who commented before, I know have a better understanding on how this all works. Also, it put me on a path to find an example of how this type of experiment is done in a clinical laboratory. You can find the study [here][1].

[1]: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4097710/#!po=50.0000

Answer: 1) From WGS data, chromosome vs episome can be distinguished by copy number (reflected in differences in read depth). With appropriate data, you can also assemble the genomes and distinguish the two by contigs.

2) Search for 'amplicon sequencing'. Note that the MiSeq, at 10M+ reads/run, is overkill unless you're barcoding 1000s of samples.
Question: I would like to extract sequences consecutively from a large fasta files and save in new fasta files. For example, I have two fasta files namely species1.fasta and species2.fasta as given below,

species1.fasta
>tag1
TGCAAGAG
>amylase
ATCGAGAT

species2.fasta
>prot8
ATGTTTGCGA
>lip1
TGATATAGAT
and I need to extract first positioned sequences from both `species1.fasta` and `species2.fasta` and save in new file namely `group1.fasta`. Same manner, the second positioned sequences to be extracted and saved in `group2.fasta`, and the `expected output` should be like this,

group1.fasta
>tag1
TGCAAGAG
>prot8
ATGTTTGCGA

group2.fasta
>amylase
ATCGAGAT
>lip1
TGATATAGAT

Like this I have to do for 3000 sequences, manual process is impossible, therefore, I have tried the following script, wherein I could not loop the process. Therefore, please help me in this regard.

my script
for file in *.fasta do cat "$file" \ | awk '/^>/ {printf("%s%s\t",(N>0?"\n":""),$0);N++;next;} \ {printf("%s",$0);} \
END {printf("\n");}' \ | awk '{if ( NR==1 ) print $0;}' \ | awk '{print $1"\n"$2}' | uniq > "new_${file}"
Thank you.

Answer: A fairly naive approach (code not tested as I'm not at a computer to do so at the moment, so it will likely be wrong in some places, but the general approach should work)

from Bio import SeqIO
import sys

for file in sys.argv[1:]:
for i, rec in enumerate(SeqIO.parse(file, 'fasta')):
file = open('./sequences_' + str(i), 'a')
file.write(rec.format('fasta'))
file.close()

This should iterate all records in all files, and write a new file with the number. If the file doesnt yet exist, it will create it, if it does, it will append to it.

Run as:

`python scriptname.py /path/to/files/*.fasta`
Question: Here, the abundance measure can represent concentration, molecular number or ppm, not in the chip-matter.

I don't know why to log them in the classification or regression model as y variable?

Could anybody help me? Thank you.

Answer: Apart from symmetry-fying the increase/decrease scale, I think taking logarithm makes the distribution gaussian(ish) for which you can do parametric tests, taking log also to some extent, takes care of heteroscedasticity (non-uniform mean dependent variance)
Question: I know for the reverse procedure, there are lots of methods available. But how can I change my fasta file with sequences like:

>gi|189094002
MDDAELNAIRQARLAELQRNAAGGGSSTNPSSSSSGGAQDSAQENMTITILNRVLTNEARERLSRQQTKITFNRKNIAADDDEDDDDFFD
>gi|68485955
MDDAELNAIRQARLAELQRNAAGGGSSTNPSSSSSGGAQDSAQENMTITILNRVLTNEARERLSRVKIVRRDSQQKQQTKITFNRKNIAGDDEDDDDFFD

to like:

>gi|189094002
MDDAELNAIRQARLAELQRNAAGGGSSTNPSSSSSGGAQDSAQENMTITILNRVLTNEARERLSR
QQTKITFNRKNIAADDDEDDDDFFD
>gi|68485955
MDDAELNAIRQARLAELQRNAAGGGSSTNPSSSSSGGAQDSAQENMTITILNRVLTNEARERLSR
VKIVRRDSQQKQQTKITFNRKNIAGDDEDDDDFFD

I know its just about inserting line breaks in each sequence, but is there any way that I can check whether the fasta file is single line fasta or not, and if is, add line breaks in each sequence..

Ps. I am using windows..

Thanks for your consideration

Answer: if your data can fit in memory, here is a simple HTML form that should do the job:

https://gist.github.com/lindenb/d420616f302a77050121ba75737c66c6
Question: Hi,
I have downloaded sra file from NCBI, SRR4242282.sra and I got three fastq files after use fastq-dump to extract fastq files from sra files.
command :fastq-dump --split-3 --gzip SRR4242282.sra
I have no idea with this result, I haven't met this before. Any suggestion would be appreciated.

Answer: Please see my previous [post][1]: this is due to the fact that BOTH paired and unpaired reads are included in these sra files. Using the `--split-files` option does not work since this would lead to fastq-files that are incomplete. What you did is correct; simply use the files ending with `_1` and `_2` will do. The remaining files contains the unpaired reads, and can be trashed.

[1]: https://www.biostars.org/p/185429/
Question: Dear all,

I am working on an organism that has 10 copies of rRNA genes on its genome. I would like to see the percentage of reads that are aligned to ribosomal RNA.

After alignment, rRNA reads align to every and each of these genes and if I perform the counting for reads that are not uniquely aligned to the genome, each one of these reads are counted 10 times. How can I count each read only once?

Any idea is highly appreciated.

Thanks

Answer: Hi, I am not sure whether this is the option that works best for you. But, I suppose you can get what you want using [Tagdust2][1]. Basically you provide the sequences of rRNA genes (I usually generate this using [BioMart][2]). I have used this in the past to remove rRNA reads before genome alignment. Run Tagdust2 and then take a look at the TagDust2 output summary file. If you do this following math you will get the number you want.

rRNA reads = total input reads - successfully extracted

When you do this, these rRNA reads will be saved on to a separate .fq file which you can use for any other downstream analyses you want to perform. May be an advantage for you.

To run Tagdust2:

> tagdust -t 8 -1 R:N -o /path/to/output/ -ref /path/to/rRNA.txt -fe 2 /path/to/input.fastq

I hope this will help!

[1]: http://tagdust.sourceforge.net/
[2]: https://m.ensembl.org/biomart/martview/8ede409a16af70d56987a76ec8ba2422
Question: I find that in hg19 fasta file, the sequence contains capital letters(ATGC) and lowercase letters(atgc), so I am wondering is there any difference between these 2 kinds of letters, or is there any special meanings?

Answer: See http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/

> Repeats from RepeatMasker and Tandem Repeats Finder (with period of 12 or less) are shown in lower case; non-repeating sequence is shown in upper case.