cannin/biostars_qa huggingface dataset:
Question: Hi all,

I am using tess3r package to plot bar plots of the Q matrix of ancestry coefficients that I got from sNMF analysis.

But I am interested to order the matrix based on counties first and then based on Q, in another word, I want to have the matrix ordered by countires and within the conuty ordered by Q.

For example this is a q matrix, I've added the county for each individuals

1.00E-04	0.75453	1.00E-04	0.24507	1.00E-04	1.00E-04	CM04487	Benin
1.00E-04	0.880116	0.0482435	0.0713402	1.00E-04	1.00E-04	CM04489	Mali
0.0822559	0.805366	0.00713594	1.00E-04	0.0198168	0.0853258	CM04493	Benin
0.967625	1.00E-04	1.00E-04	1.00E-04	1.00E-04	0.0319755	CM06496	Burkina Faso
0.83507	0.0272506	1.00E-04	0.00335877	1.00E-04	0.134121	CM06501	Burkina Faso
0.435478	0.235613	1.00E-04	0.131601	0.0401917	0.157016	CM06505	Guinea
0.302861	0.313093	0.0068676	1.00E-04	0.00161187	0.375466	CM06510	Mali
0.335162	0.365021	0.0119398	0.0315991	0.0107296	0.245548	CM06511	Ghana
0.502391	0.0369072	1.00E-04	0.0648057	0.0577006	0.338095	CM06513	Burkina Faso
1.00E-04	0.734269	0.000945988	0.259324	0.00526157	1.00E-04	CM07885	Guinea

in tess3r package running this part won't give what I want, because it sorts based on Q only.

barplot(Q.matrix, sort.by.Q = FALSE,
border = NA, space = 0,
col.palette = my.palette,
xlab = "Individuals", ylab = "Ancestry coefficients") -> bp

Guid and help is much appriciated.

Thanks!

Answer: Order by column 8 first and then "sub"-order by column 7

Q.matrix <- Q.matrix[order(Q.matrix[,8],Q.matrix[,7]),]

EDIT: Actually this might be wrong... what column is Q? If it's not a variable in your data.frame it will be hard to use it for ordering.

Otherwise you can split your data.frame based on country with `split(Q.matrix,Q.matrix[,8])` and then plot each country separately (using `lapply` or a `for` loop).
Question: Is there a way to retrieve an article's citation index via Pubmed or Pubmed's eutils or any other way in Pubmed/Entrez? Or a way to filter the search results based on citation index in pubmed?

Answer: If by "citation index" you mean the number of times a paper is cited, then not really. The closest you'll come is the number of citations indexed in Pubmed Central. Using R:

library(rentrez)

get_pmc_cites <- function(pmid){
res <- entrez_link(db="pmc", dbfrom="pubmed", id=pmid)
return(length(res$pubmed_pmc_refs))
}
get_pmc_cites(20203609)
# [1] 25

***EDIT***

Just noticed the summary xmls for pubmed have the PMC reference count, which might make multiple calls quicker:

recs <- entrez_summary(db="pubmed", id=c(20203609, 20203610))
sapply(rec, "[[", "PmcRefCount")
# [1] 25 200
Question: Hi, I was trying to use vg gamsort to sort and index the gam file so as to use vg chunk later. I get the following error when I run gamsort

[E::bgzf_flush] File write failed (wrong size)

[E::bgzf_flush] File write failed (wrong size)

terminate called without an active exception

I think this is because it is running out of space in the temporary directory. I am not sure where this tmp directory is or if I have access to it, as I work in a shared cluster. Is there a way to specific the tmp directory as it is possible for vg index? I got the same error while running vg index but was able to rectify it using the -b option. While vg index allows indexing of a sorted gam file, mine comes directly from vg map so don't think it is sorted.

Thanks!

Answer: `vg gamsort`, as well as most other vg commands, takes the temporary directory from environment variables. The primary variable is `TMPDIR`, but vg also checks a number of other variables (`TMP`, `TEMP`, `TEMPDIR`, `USERPROFILE`) if that has not been set. If everything else fails, vg falls back to `/tmp`.

The main exception are temporary files created by external libraries. We try to make the libraries use the same directory vg uses, but that's not always possible.
Question: I am currently using some genes in Affymetrix format, like:
- 200007at
- 200011sat
- 200012x_at

What do "at", "sat", and "x_at" stand for?

Thanx

Answer: Here is an answer from the Affymetrix FAQ:

The last three characters identify the probe set strand. Probe sets that are designed to detect the anti-sense strand of the gene of interest are annotated with "_at". Probes in an identical probe set (_s set) all cross-hybridize to the same set of sequences that are used for the design (including any additional pruning sequences if provided). Probes in a mixed probe set (_x set) contain at least one probe that cross-hybridizes with other sequence(s) used for the design.

See tjhe above link for more information, including a nice figure depicting these relationships and definitions.
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: Some online resources:

Quick-R
a set of tutorials that I found to be to be a very lucid and practical introduction for general data analysis with R.

Data Analysis Examples at UCLA
a nice matrix of various canonical analytical tasks in R as well as SAS, SPSS, etc.

R and BioConductor Manual
a concise introduction to R for bioinformatics specifically
Question: The userguide says:
![userguide][1]

But I find another explanation saying that，the -Z should be set to the size of genome times 2.

I am in a mess now. How should I set the -Z , and I also want to know what's the E-value?

[1]: /media/images/092d9bed-d881-41c3-9105-d1c0639c

Answer: To calculate Z please visit this [link](https://docs.rfam.org/en/latest/genome-annotation.html) and use the function

```
esl-seqstat my_reference.fasta*
```
* reference fasta should be your genome.

E-score:
The E-value is the statistical significance of the hit: the number of hits we’d expect to score this highly in a database of this size (measured by the total number of nucleotides) if the database contained only nonhomologous random sequences. The lower the E-value, the more significant the hit.

Z: This option ensures that the reported E-values are accurate.

I hope that his helps:
Question: Hi all,

Does somebody know which package can be used for plots like this?

> The dot plots illustrates the number of detected fusion transcripts per megabase per sample normalized by the sequencing coverage. Tumor types were sorted according to the fraction of samples with fusions.

I am trying to recreate this plot from the gene fusions paper.
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4468049/

The data is present here:
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4468049/bin/NIHMS632238-supplement-2.xlsx
(Table S1) and using only columns "Cancer, FileName,FusionPair"

I have tried a couple of options from ggplot package but these are different plots.
http://www.sthda.com/english/wiki/ggplot2-dot-plot-quick-start-guide-r-software-and-data-visualization

Thanks

Ron

Answer: Here is the start:

library(ggplot2)

#dummy data
mydat <- mtcars[ , c("cyl", "mpg")]
mydat$Cancer <- paste0("cancer_", mydat$cyl)

#Order by cancer and mpg (number of fusions per megabase)
mydat <- mydat[ order(mydat$Cancer, mydat$mpg), ]
mydat$x <- seq(nrow(mydat))

# make pretty x labels one per cancer
mydat$xlabel <- c(1, diff(as.numeric(as.factor(mydat$Cancer))))
mydat$xlabel <- ifelse(mydat$xlabel == 1, mydat$Cancer, NA)

# plot
ggplot(mydat, aes(x, mpg, col = Cancer)) +
geom_point() +
scale_x_continuous(breaks = which(!is.na(mydat$xlabel)),
labels = mydat$xlabel[ !is.na(mydat$xlabel)])
Question: Hi!

I want to reverse translate a short protein sequence (around 50 aa) to DNA sequence using all possible codons for all the positions. If Glu is GAA and GAG, I want to have a sequence variants with both.
So as a result I want to have a collection of all possible DNA sequences that could have made this protein (regardless the codon usage).

Is there any tool like this available? If not what's the fastest way to obtain such file?

Thank you!

Answer: Interesting question - I never thought of it that way. I ended up with something like this:

import itertools

d = {
'A': ['GCA', 'GCC', 'GCG', 'GCT'],
'C': ['TGC', 'TGT'],
'D': ['GAC', 'GAT'],
'E': ['GAA', 'GAG'],
'F': ['TTC', 'TTT'],
'G': ['GGA', 'GGC', 'GGG', 'GGT'],
'H': ['CAC', 'CAT'],
'I': ['ATA', 'ATC', 'ATT'],
'K': ['AAA', 'AAG'],
'L': ['CTA', 'CTC', 'CTG', 'CTT', 'TTA', 'TTG'],
'M': ['ATG'],
'N': ['AAC', 'AAT'],
'P': ['CCA', 'CCC', 'CCG', 'CCT'],
'Q': ['CAA', 'CAG'],
'R': ['AGA', 'AGG', 'CGA', 'CGC', 'CGG', 'CGT'],
'S': ['AGC', 'AGT', 'TCA', 'TCC', 'TCG', 'TCT'],
'T': ['ACA', 'ACC', 'ACG', 'ACT'],
'V': ['GTA', 'GTC', 'GTG', 'GTT'],
'W': ['TGG'],
'Y': ['TAC', 'TAT'],
'_': ['TAA', 'TAG', 'TGA'],
}

def generator(protein):
l = [d[aa] for aa in protein]
for comb in itertools.product(*l):
yield "".join(comb)

if __name__ == '__main__':
import sys
protein_seq = sys.argv[1]
g = generator(protein_seq)
for dna_seq in g:
print(dna_seq)

Run:

python script.py MKS

Output:

ATGAAAAGC
ATGAAAAGT
ATGAAATCA
ATGAAATCC
ATGAAATCG
ATGAAATCT
ATGAAGAGC
ATGAAGAGT
ATGAAGTCA
ATGAAGTCC
ATGAAGTCG
ATGAAGTCT
Question: I have a bunch of human Ensembl transcript ids from an old release, that I want to find the ids in the current version.

For example, the transcript: [**ENST00000360913**][1], is was removed since release 76. But If I BLAST its sequence, I can find a 100% match for another transcript in the current release (80): **ENST00000585580**.

Is there a way to do it automatically?

Is there a resource of mapping such ids?

[1]: http://www.ensembl.org/Homo_sapiens/Transcript/Idhistory?t=ENST00000360913

Answer: The two transcripts are distinct transcripts that code for the same protein (see the archived [Ensembl browser 75][1]), meaning what you are doing isn't strictly mapping it to a newer version of the transcript, but to an entirely different transcript.

If you want to do this programmatically, you could probably implement a solution using [biomaRt][2] in R. I've put together a solution that works for the example you've given, though I haven't tested it with any others. It should give you an idea of how to go about this, if nothing else:

```r
library(biomaRt)

#transcript IDs to "convert"
transcript.ids <- "ENST00000360913"

# define biomart objects, current and archived
mart.current <- useMart(biomart="ensembl",dataset="hsapiens_gene_ensembl")
mart.archive <- useMart(host="feb2014.archive.ensembl.org",biomart="ENSEMBL_MART_ENSEMBL",
dataset="hsapiens_gene_ensembl")

# query biomart for the CCDS ID, shared by the two transcripts
archive.results <- getBM(attributes=c("ccds","transcript_biotype","ensembl_transcript_id"),
filters="ensembl_transcript_id",values=transcript.ids,mart=mart.archive)
current.results <- getBM(attributes=c("ensembl_transcript_id","transcript_biotype","ccds"),
filters="ccds",values=archive.results$ccds,mart=mart.current)

mapped <- merge(archive.results,current.results,by=c("ccds","transcript_biotype"),
suffixes=c("_archived","_current"))```
```

[1]: http://feb2014.archive.ensembl.org/Homo_sapiens/Transcript/Sequence_Protein?db=core;g=ENSG00000178093;r=19:19623230-19626838;t=ENST00000360913
[2]: http://bioconductor.org/packages/release/bioc/html/biomaRt.html
Question: Hi All,
I need some help interpreting my DEseq2 results. Here is a bit of my DEseq2 code:

dds$cond <- relevel(dds$cond, ref = "C")
dds = DESeq(dds, fitType = "parametric")
my.names<-resultsNames(dds)
my.names
[1] "Intercept" "cond_T_vs_C"

Then I did:

res<- results(dds, contrast=list( "cond_T_vs_C"), test="Wald")

Now looking at the result, for this particular gene I have:

baseMean log2FoldChange lfcSE stat pvalue padj
Seq_483975_2034 1995.245 -1.354257 0.4045764 -3.347346 0.000815893 0.005723768

Can someone please explain me whether this gene is down regulated in Control (C) or in Treatment (T)? Thanks

Answer: This means treatment is lesser than control by a log fold change of 1.354. You can check this by

counts(dds, normalized=TRUE)

you can see the normalized counts of control and treatment.

Aifu.
Question: Hi all,

I have done BLASTp searches on multiple genomes and pfam analysis. So now, I have numerous .pfam result files (text files) and I would like to extract only identical matches between two genomes.

Instead of doing it manually, is there a bash command or else to **extract all identical words from text files**?

I tried this:

comm -12 file1.txt file2.txt > output.txt

But as I expected, I only have the common lines between the files while I need all common words...

If you have any ideas, please let me know!

Thank you in advance for your help.

***EDIT: file examples***

**Here is an example of file1:**

# --- full sequence ---- --- best 1 domain ---- --- domain number estimation ----
# 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
#------------------- ---------- -------------------- ---------- --------- ------ ----- --------- ------ ----- --- --- --- --- --- --- --- --- ---------------------
WxL PF13731.7 Lactococcus_893 - 3.7e-07 30.4 7.2 5.2e-07 29.9 7.2 1.2 1 0 0 1 1 1 1 WxL domain surface cell wall-binding
DUF916 PF06030.13 Lactococcus_894 - 7.3e-20 71.3 1.3 8.7e-20 71.0 1.3 1.0 1 0 0 1 1 1 1 Bacterial protein of unknown function (DUF916)
DUF3324 PF11797.9 Lactococcus_895 - 3.2e-31 108.2 1.3 3.9e-31 107.9 1.3 1.1 1 0 0 1 1 1 1 Protein of unknown function C-terminal (DUF3324)
APH PF01636.24 Lactococcus_896 - 3.2e-07 30.5 0.1 4.4e-07 30.0 0.1 1.1 1 0 0 1 1 1 1 Phosphotransferase enzyme family

**Here is an example of file 2:**

# --- full sequence ---- --- best 1 domain ---- --- domain number estimation ----
# 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
#------------------- ---------- -------------------- ---------- --------- ------ ----- --------- ------ ----- --- --- --- --- --- --- --- --- ---------------------
Beta-lactamase2 PF13354.7 Lactococcus_549 - 1.9e-33 115.8 0.0 2.8e-33 115.2 0.0 1.3 1 0 0 1 1 1 1 Beta-lactamase enzyme family
S1 PF00575.24 Lactococcus_550 - 0.014 15.8 0.7 0.47 10.9 0.7 2.3 1 1 0 1 1 1 0 S1 RNA binding domain
DUF3324 PF11797.9 Lactococcus_551 - 3.4e-08 33.2 0.5 3.4e-08 33.2 0.5 1.9 2 0 0 2 2 2 1 Protein of unknown function C-terminal (DUF3324)

As in the two .txt files there is only the DUF3324 target name in common, I would want in an output file simply like this:

DUF3324

I'm really thankful if you want to help me on this...

Answer: The command you found is correct, you just need to use it correctly:

comm -12 <(tail -n+4 file1.txt | cut -d ' ' -f 1| sort ) <(tail -n+4 file2.txt | cut -d ' ' -f 1| sort )

You don't need to do all of this in one line if you don't want to, but you do need to:

1. Trim the header lines out (`tail -n+4`)
2. Isolate just the first column, as defined by whitespace (`cut -d ' ' -f 1`)
3. `sort` the files in order for `comm` to work (if they aren't already sorted, which these aren't).
Question: Hi There,

I did copynumber call and performed segmentation using DNAcopy.

```r
cnP1 = read.table ("copyNumber_varscanIdhP1.copynumber", header=TRUE)
CNA.object <- CNA( genomdat = cnP1[,7], chrom = cnP1[,1], maploc = cnP1 [,2], data.type = 'logratio')

CNA.object.smoothed <- smooth.CNA(CNA.object)
CNA.object.smoothed.seg <- segment(CNA.object.smoothed, verbose=0, min.width=2)
segmented = CNA.object.smoothed.seg$output

write.table (segmented[,2:6], file="cbsP1", row.names=F, col.names=F, quote=F, sep="\t")
```

Now, I want to merge these output segments:

```
perl mergeSegments.pl
USAGE: mergeSegments.pl [segments] OPTIONS
segments - A segments file with p-values from the DNAcopy library. This should be tab- or space-delimited
with a header and the following columns: ID, chrom, loc.start, loc.end, num.mark, seg.mean, bstat, pval, lcl, ucl.
OPTIONS:
```

WHAT is the input file here?

Is the input to the script is output from DNAcopy output.

But the columns are different:

```
head cbs_P1
chr1 567331 569882 6 -2.3721
chr1 761847 905536 162 -0.1232
chr1 905636 3549990 2846 -0.0335
```

Can you please suggest me to make the input file in compatible format.

Thanks!

**Edit**: This is how you could get input file (output from DNACopy) to mergeSegments.pl

```r
seg.pvalue <- segments.p(CNA.object.smoothed.seg, ngrid=100, tol=1e-6, alpha=0.05, search.range=100, nperm=1000)
write.table (seg.pvalue, file="blah", row.names=F, col.names=F, quote=F, sep="\t")
```

Answer: Its the output from segments.p function from DNACopy.
Question: Hi all,

I have a question regarding coding regions in bacterial genomes: **can genes be encoded on both strands like this ?**

![enter image description here][1]

[1]: /media/images/fceb469c-9a40-47f9-969b-fdce9ddb

**If so, which one is considered as the sense strand ?**

Thank you so much for your help!

Answer: yes, for example ryeA/ryeB on E Coli K12:

[link to NCBI view][1]

[1]: https://www.ncbi.nlm.nih.gov/nuccore/48994873?report=graph&tracks=[key:sequence_track,name:Sequence,display_name:Sequence,id:STD649220238,annots:Sequence,ShowLabel:false,ColorGaps:false,shown:true,order:1][key:gene_model_track,name:Genes,display_name:Genes,id:STD3194982005,annots:Unnamed,Options:MergeAll,CDSProductFeats:false,NtRuler:true,AaRuler:true,HighlightMode:2,ShowLabel:true,shown:true,order:4][key:feature_track,name:Repeat%20region,display_name:Repeat%20region,id:STD3463812800,subkey:repeat_region,annots:Unnamed,shown:true,order:5]&assm_context=GCA_000005845.1&v=1920965:1921462&c=C0C0C0&select=null&slim=0

> If so, which one is considered as the sense strand ?

I think it's meaningless
Question: Dear all,

to test the alignment procedure, I inserted a genome in reverse sense into an artificial set of fastq files and aligned against the reference genome in normal orientation. The resulting alignment does not show any mapping.

Is it normal that the aligner could not map the reads because they were in the reverse orientation? I thought that orientation could not affect the alignment procedure...

Thank you

Answer: Reference genome orientation of course matters. What is interchangable would be the reference genome in its standard representation and its reverse complement but not its reverse. That would completely ignore that DNA has a 5'->3' polarity that cannot simply be changed at will.
Question: Is there any possibility to predict the specific transcript (obtained by alternative splicing) of a gene present in a particular individual just by exome sequencing and calling variants?

Because we got an analysis report from a recognized institute for an exome sequence we sent which says that they have identified a nucleotide substitution of xxxxx gene, A>G and they have given the transcript ID [ENST00000xxx] as well.

And when I explored further I found this from 1000genomes where they only did exome sequencing but have annotated SNPs with transcript IDs.

http://pilotbrowser.1000genomes.org/...0000323668#ALL

Please could somebody explain this to me? I'm a novice bioinformatician

Answer: On the example from 1000 genomes, if you click on one of the "Show" links, e.g. for "Frameshift coding variants", you'll see exactly how this works. Basically, a single variant can end up getting associated with all of the transcripts that overlap it. So if you have a variant at 5:149720925, then it'll have 4 associated transcripts. These sorts of annotations are mostly for convenience, so you can quickly assess how causative a variant is (if none of the associated transcripts are expressed in your tissue of interest, then the variant is unlikely to be causative).
Question: Hi everybody !

So there is my problem : I am searching for the average distance between a known motif inside sequence, and extend this to a list of sequences...
The first part is done, the second part (extend to a list of sequences) is the problematic one ! So, here the way i am doing the first part :

> source("motifOccurrence.R") #https://www.r-bloggers.com/calculate-the-average-distance-between-a-given-dna-motif-within-dna-sequences-in-r/

> library("seqinr")

> df <- readDNAStringSet("X.fasta")

> df2 <- df[[1]]

> motif <- c("T", "C", "C", "A")

> coord <- coordMotif(df2, motif)

> motidist <- computeDistance(coord)

> motidist

[1] 152

It's appear that the first sequence of my fasta list have an average distance of 152 nucelotide between two TCCA motif.
And, i don't know how automatize this to all my list in df...

Thanks by advance for the help.

Kévin

Answer: The [Rho][1] function in the seqinr package calculate the indipendence of occurence of a k-words using the expected word count.
We can modify it to be run just for our target:

library("seqinr")
data(ec999)

exp_count<-function (sequence, wordsize = 2, alphabet = s2c("acgt"))
{
wordcount <- count(sequence, wordsize, freq = FALSE, alphabet = alphabet)
uni <- count(sequence, 1, freq = TRUE, alphabet = alphabet)
expected_wordfreq <- function(wordsize, uni) {
if (wordsize == 1)
return(uni)
else kronecker(uni, expected_wordfreq(wordsize - 1, uni))
}
expected_wordcount <- sum(wordcount) * expected_wordfreq(wordsize,
uni)

names(expected_wordcount)<-words(wordsize,alphabet = alphabet)
return(expected_wordcount)
}
#exp_count(ec999[[1]],3)

[1]: https://rdrr.io/rforge/seqinr/man/zscore.html
Question: I am using rsem to calculate expression levels in a dataset of paired end reads. I want to use STAR as the aligner. I have written the command like this:

rsem-calculate-expression --star -p 16 --paired-end DT_01_R1.fastq DT_01_R2.fastq /path/to/reference/genome/ rsem_output/DT_01

I'm pretty sure this will work, What I don't know is this: do I need to generate an annotated reference genome for STAR before running this? I'm using the zebrafish genome, and I have downloaded both the .fasta file and the .gtf file from Ensembl. What do I need to do with these before I run the aligning and expression calculation I've written above? Or do I just point it to the directory containing the .fasta and .gtf files? Thanks very much!

Answer: You have to prepare the transcriptome reference first, see [rsem-prepare-reference][1]. See a tutorial [here][2].

Alternatively, you can map with STAR to the genome (you will have to prepare the index yourself), using `STAR --quantMode transcriptomeSAM`, and then use STAR's Aligned.toTranscriptome.out.bam output as input to `rsem-calculate-expression`.

[1]: http://deweylab.biostat.wisc.edu/rsem/rsem-prepare-reference.html
[2]: https://github.com/bli25broad/RSEM_tutorial
Question: Hi all!

Are the values in geo .soft files log2 transformed? are they normalized? Thanks in advance

Answer: Yes, some of them are log2 transformed, some are normalized, some are neither. Presuming you want to know which, the only way to know with certainty is to read the meta information at the top of the .soft file.