# Question
Title: Is the Bowtie2 MAPQ score dependent on the database index?
I have two very similar genomes in my bowtie2 database (nuc similarity \> 90 %). When mapping with bowtie2 against this database, most of the reads that are mapped to one of these genomes have a MAPQ score of 0. Only at the positions that are different between the two genomes I get higher MAPQ scores.

I also created a database with just one of the genomes, inspecting the MAPQ score of a single read. This single read gets a MAPQ score of 24 when the database contains only one genome, but drops to zero when two similar genomes are in the database.

Is this the expected behavior? My understanding of the MAPQ score computation is, that the alignment score of a read at a specific location is only compared to the score at other positions within the same reference, not across the whole database. If the alignment scores for a read are then pretty close to each other, it is assumed that the read does not map uniquely, which is reflected in a low MAPQ score.

So my question is, whether the alignment score of a read is compared to other alignments of this read within the same reference OR against alignments with the whole database index? Changing the database, would then result in different MAPQ scores.

# Answer
Yes, this is expected. I am not sure what you mean by "reference" vs "database", but what you do is to index a fasta file or a collection of fasta files. That index you create by that is what `bowtie2` alignes against. It does not know whether the index consisted of one, ten or a thousand fasta files, it just maps against the entire index and then check whether it can uniquely (or not) place the reads. If you combined two genomes into one index and the genomes are very similar then indeed many reads will map equally well to them and that manifests as MAPQ 0 (or low).
> 1 votes

# Answer
Yes, this is expected.

If a read maps equally well to two different locations in the reference database, then the minimum probability that read maps to another location is 0.5. MAPQ scores are ten times the negative `log10` of the probability of a mismatch, so a mismatch probably of 0.5 would lead to a *maximum* MAPQ score of $-10\*log\_{10}(0.5) = 3$. If there's any doubt about the mapping (e.g. a mismatch in the sequenced read), or the read could map to more than two locations, then that probability will drop further. I would not be surprised if these MAPQ scores, where the location is impossible to precisely determine, were rounded down to 0.
> 0 votes

---
Tags: bowtie2, mapq
---


# Question
Title: Differentially expressed genes analysis in Seurat
For the differentially expressed genes analysis, is it possible to check for DEGs based on the levels already identified in the object? For example, my dataset contains cells from 11 subjects, which are marked as so (orig.ident). I want to analyze the differentially expressed markers between them. Will it be accurate to do so? Or should I follow the clustering analysis in Seurat? I want to see the difference in genes between the cells according to the subjects but I'm afraid that may not be the right way to do so.

```
  head(BC@meta.data)
         orig.ident nCount_RNA nFeature_RNA percent.mito
 BC01_02       BC01   998944.0         9189 0.2387108038
 BC01_03       BC01   999696.4         9548 0.2925345423
 BC01_04       BC01   999057.5         7440 0.0009893724

```

Is clustering analysis absolutely necessary to see the DEGs between them? or can I directly move to FindMarkers based on the levels (orig.ident)?

# Answer
> 1 votes
@swbarnes2 code is one correct way to find marker genes that differ between samples. You could also change the `Idents` slot and then just use `FindAllMarkers`.

It's not clear what you mean by accurate in your comment. If you look for marker genes between samples (orig.ident) without clustering, Seurat will use expression data from all the cells attributed to each sample to find sample-specific markers. This approach can be right or wrong depending on the question you are asking.

For example, this paper asked which genes were unique to each patient's tumor (see Figure 2D). To recreate their analysis, you would restrict your Seurat object to only include tumor cells (removing other cell types like immune cells and fibroblasts) and then perform `FindMarkers` on sample origin. This answers which genes are specifically expressed on each patient's tumor cells, averaged over the different tumor cell subpopulations (in HNSC there are basal, differentiated, stressed, and pEMT subpopulations). This analysis could be a good idea if you are trying to find a sample-specific marker that is expressed on the majority of tumor cells in a sample. The downside to this approach is that it averages across tumor cell subpopulations. If there is a rare tumor subpopulation (say a set of cancer stem cells) that do not express the same markers as the majority of cells, the markers Seurat finds may not be expressed on this rare population.

It doesn't make sense to look for sample specific markers when your Seurat object includes several different types of cells (epithelial cells, immune cells, stromal cells etc). In other words, you should probably restrict your analysis to some major cell type before looking for marker genes. If you don't, your results will be analogous to a bulk RNA-Seq differential expression, and will be significantly influenced by the cell type composition between samples which defeats the purpose of doing single cell RNA-sequencing. In this case, one could consider the results inaccurate, as it will be difficult to attribute marker genes to one cell type.

Another way to interpret your question about accuracy is whether calling marker genes between samples provides accurate statistical inference (are the marker genes identified truly up/down in specific patients). If you are using single cell data to compare expression between conditions (normal vs tumor, control vs treated cells), it is advised to not use Seurat FindMarkers or similar tools when you have multiple samples. The problem is that FindMarkers treats each cell like its own sample, even though all cells from a sample are similar and thus correlated. This violates the IID assumption used by the Wilcox test and other methods used by FindMarkers and results in really low p-values and inflated false positive rates. In this scenario, recent studies suggest you should use a pseudo-bulk approach to prevent a high number of false positives. muscat is one method that performs pseudo-bulk analysis, and its paper explains pseudo-bulk analysis in more detail.

**TL;DR** If you are going to look for sample-specific (or other metadata-based) marker genes, it's probably a good idea to first restrict the analysis to a major cell type. If you are using single cell data to perform differential expression analysis comparing conditions, use pseudo-bulk approaches.

# Answer
> 0 votes
I'm not quite sure what you are asking, but you can use anything in the metadata as a grouping with FindMarkers.

```
DefaultAssay(combined.all) <- "RNA"
Idents(combined.all) <- "orig.ident"

for (sample in unique(Idents(combined.all))) {
  clustermarkers <- FindMarkers(combined.all, ident.1 = sample)
  print(paste("Cluster", sample))
  print(head(clustermarkers, n = 10))
}

```

---
Tags: r, scrnaseq, seurat, differential-expression
---


title: Imputing small region of the genome

text:
# Question
Title: Imputing small region of the genome
If I'm looking for a specific SNP in my SNP-Chip data and it isn't there, are there any tools that let me quickly impute that SNP from surrounding SNPs rather than running a lengthy 'whole chromosome' imputation job?

If so, roughly how many upstream and downstream SNPs are generally required to predict the haploblock 'accurately'?

# Answer
> 3 votes
I would take something like 1.5Mb either side of the SNP of interest (so a 3Mb chunk), subset that out from the target and reference panel and then perform imputation on that chunk alone. I think that size will give you more than enough LD, since LD decays by about 500kB in most human populations by then.

shapeit4 has an option to just impute a single region:

```
shapeit4 \
    --input unphased.bcf \
    --map chr20.b37.gmap.gz \
    --region chr21:1231213-1413143 \ ## change this bit
    --reference reference.bcf \
    --output phased.bcf

```

---
Tags: snp, imputation, haplotypes, tools
---


# Question
Title: Select synonymous sites from a multiple sequence alignment
Could someone kindly recommend a tool or R package that can identify synonymous sites in a multiple sequence alignment?

I wish to select those taxa for tree reconstruction and other downstream analyses within subset of the alignment.

To respond to the questions below: yes, I am interested in identifying synonymous sites within a multiple sequence alignment (DNA) of ~70 virus samples with limited genetic diversity. My goal is to subset the multiple sequence alignment including synonymous sites only, because we believe that there are parallel non-synonymous mutations occurring that do not provide information on shared evolutionary history.

Thank you @M\_ for the really helpful information below. In Mega, the S tool identifies singleton sites rather than synonymous sites in a multiple sequence alignment. I'm wondering if there is another tool or a PAML option to subset the synonymous sites?

Thank you in advance and best!

# Answer
> 0 votes
I proposed Syn-SCAN (below). This was used for HIV and as you are looking at viruses seems a good fit.

**There's a simpler way** ... just make a invariable amino acid alignment, i.e. discard all amino acid variation. MEGA will do that, or use an algorithm. You then simply check the remaining nucleotide sequences for identity and discard the identity (lots of algorithms do this).

What you are left with is, is what you seek, PURE SYNONYMOUS mutations.

1. discard taxa with amino acid variation
2. retain taxa nucleotides with variation.

However, just discarding amino acid variation and retaining everything else will either by identity or synonymous. To be honest, its a virus so the amount of sequence identity between taxa is likely to be very little for protein loci.

Its really that simple. I honestly wouldn't over complicate it.

Sorry for the singleton vs. synonymous error, I haven't looked at MEGA in ages. Its a bit bizarre it discarded the only thing it was actually used for.

---

The short answer to the question is the 3rd codon position is the traditional approximation for synonymous mutation. A small proportion of 1st codon position sites are neutral and a small proportion of 3rd codon position are non-synonymous.

I agree with @terdon about and @Community regarding additional information sought, because synonymous mutations are heavily subject to compositional bias and that is a big subject. The key issue is the whether the taxa/species in the alignments are known for compositional bias and the genetic distance between different members of the alignment.

The traditional package used for synonymous mutations is/was **MEGA**, currently version 11. Albeit, recent 'upgrades' seem to remove the synonymous distance matrix calculation. HOWEVER, they may have transferred this capability onto its what appears to be upgraded GUI alignment viewer. What I think is if you:

* load the alignment into the GUI;
* click alignment (big button far left of menu);
* you will now see the alignment;
* in the middle of the menu bar are 4 buttons 'C', 'V', 'Pi' and 'S'
* I think 'S' = synonymous
* All synonymous sites in the alignment are now highlighted and there is likely a method of extracting this information.

C=constant, V=variable Pi=Nei approximation of diversity (population genetics)

If this is correct then those are the subset sought, please note however you will need to select/deselect the taxa manually and reassess the synonymous sites after each 'point and click' iteration (the S button will remain active I suspect between each selection/deselection). What you appear to want is a synonymous tree and this function I believe was removed from MEGA.

MEGA is rarely used except for teaching because its functionality is available elsewhere that permit a given algorithm to be integrated into a data pipeline, which MEGA would not permit (might have recently changed). Hence my recollection of MEGA is a little vague.

---

The package used for synonymous work is PAML and remains the undisputed standard without question in the context of formal selection analysis, which is heavily dependent on a estimate of synonymous mutations. What you refer to as 'synonymous mutation' isn't the same as the phylogenetic definition, i.e. your definition depends on the other members of the alignment, whereas this isn't an issue in the phylogenetic definition (because the tree circumvents this, i.e. delineates subsets) except if compositional bias occurs ... huge amounts of work is/has been dedicated to compositional bias.

Final point is the widely used phylogenetic library in R called `ape` will not perform a synonymous mutation calculation - I think, however phangorn library in R might do, but is more likely within the context of selection analysis. Thus to answer the question directly, I don't think R will do this calculation.

---

It transpires S= singleton (sorry) ....

Have a look at Syn-SCAN. The link is here https://hivdb.stanford.edu/pages/synscan.html . I have never used it and the download link is on the website. The good thing is its written in Perl so it would be easy to modify and check (if you can read/write Perl).

PAML: you can't detach the synonymous mutation calculation from the non-synonymous mutation calculation, without delving into the C-code (which I wouldn't recommend).

---
Tags: sequence-alignment, phylogenetics, phylogeny, synonymous-mutations
---


title: baseMean threshold

text:
# Question
Title: baseMean threshold
I have an RNA-seq dataset and I am using `DEseq2` to find differentially expressed genes between the two groups. I used pre-filtering to remove any genes that have no counts or only one count across the samples. furthermore, I want to:

* remove genes in low counts by using a base mean threshold.
* remove genes that have low counts compared to the rest of the genes.

**Question**

Is there a common threshold used for the baseMean or a way to work out what this threshold should be?

# Answer
Adding my answer from the crosspost:

I would not use the baseMean for any filtering as it is (at least to me) hard to deconvolute. You do not know why the baseMean is low, either because there is no difference between groups and the gene is just lowly-expressed (and/or short), or it is moderately expressed in one but off in the other. The baseMean could be the same in these two scenarios. If you filter I would do it on the counts. So you could say that all or a fraction of samples of at least one group must have 10 or more counts. That will ensure that you remove genes that have many low counts or zeros across the groups rather than nested by group, the latter would be a good DE candidate so it should not be removed. Or you do that automated, e.g. using the edgeR function filterByExpr.
> 1 votes

---
Tags: r, rna-seq, deseq2, gene-expression
---


# Question
Title: Pipeline for paired end RNA sequence data to proteins
Forgive the basic question here, but I'm super novice at this ...

I have a series of paired-end RNA fastq files (e.g. `SRR10720226_1.fastq.gz` and `SRR10720226_2.fastq.gz`). I'm trying to figure out what software I need to use to translate those into proteins.

I was trying to use `biopython` to do this using the `translate` method. The problem is `biopython` doesn't seem to consider both files together, only one `fastq` at a time.

Which tool should I use to work with the reads together? Is there anything else I should be doing besides just translating to protein?

# Answer
> 1 votes
This is a much more complicated question than it might seem.

First, you need to understand how RNA-Seq works, and what your data really is. Your "paired-end files" contain reads, which will contain fragments of transcripts.

Since those are only fragments of transcripts, you don't know what frame they are in, so you don't have a single possible translation. And each read would only give you a small part of the protein.

Since the data is for humans, your best course of action would be to first, align the reads to the annotated human genome (or transcriptome). Then, count the number of reads for each gene (or transcript), and decide on a condition to call a gene (or transcript) present in the sample (typically by enrichment between two sets of samples, in your case a threshold might be good enough). Finally, for the genes (or transcripts) present in the sample, find out what protein (or isoform) they correspond to, and if needed retrieve the sequence.

To be clear, you will have to learn a lot more to get that done. For the first step, easier tools might include kallisto or salmon. You would probably want to use biological replicates; SRP237926 seems to have several samples, reading the paper would be necessary to determine whether they are useful.

---
Tags: python, fastq, rna
---


title: Error in DESeq2 analysis with dimnames

text:
# Question
Title: Error in DESeq2 analysis with dimnames
I'm getting the following error in a DESeq2 analysis and I can't seem to figure out the issue:

`Error in `assays<-`(`*tmp*`, withDimnames = withDimnames, ..., value = `*vtmp*`) : please use 'assay(x, withDimnames=FALSE)) <- value' or 'assays(x, withDimnames=FALSE)) <- value' when the dimnames on the supplied assay(s) are not identical to the dimnames on RangedSummarizedExperiment object 'x'`

Here is my code:

```
library(DESeq2)
library(edgeR)

url <- "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file"
utils::download.file(url, destfile="GSE63310_RAW.tar", mode="wb")
utils::untar("GSE63310_RAW.tar", exdir = ".")
files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt", "GSM1545538_purep53.txt",
           "GSM1545539_JMS8-2.txt", "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt",
           "GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt", "GSM1545545_JMS9-P8c.txt")
for(i in paste(files, ".gz", sep=""))
  R.utils::gunzip(i, overwrite=TRUE)

files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt",
           "GSM1545538_purep53.txt", "GSM1545539_JMS8-2.txt",
           "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt",
           "GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt",
           "GSM1545545_JMS9-P8c.txt")

x <- readDGE(files, columns=c(1,3))

samplenames <- substring(colnames(x), 12, nchar(colnames(x)))

colnames(x) <- samplenames
group <- as.factor(c("LP", "ML", "Basal", "Basal", "ML", "LP",
                     "Basal", "ML", "LP"))
x$samples$group <- group

#DESeq2 Analysis
phenodata <- new("AnnotatedDataFrame", data = data.frame(celltype = x$samples$group))
counts <- x$counts
ids <- c(1,3,4,6,7,9)

tmpcounts = counts[, ids]
tmpph = as(phenodata[ids, ], "data.frame")
rownames(tmpph) <- colnames(tmpcounts)

dds = DESeqDataSetFromMatrix(countData = tmpcounts, colData = tmpph, design = ~ tmpph$celltype)

```

Here are my package versions

```
> packageVersion("DESeq2")
[1] ‘1.30.1’

> packageVersion("SummarizedExperiment")
[1] ‘1.20.0’

```

# Answer
What works is:

```
dds = DESeqDataSetFromMatrix(countData = data.frame(tmpcounts), colData = tmpph, design = ~celltype)

```

Mind the `data.frame()` around the count data object (which was a matrix) and removal of `tmpph$` for the design because design should reference one or more columns of the colData. In newer versions of DESeq2 that by the way works without the `data.frame()` but I am not sure why. You might want to dig that out yourself, looking at the source code of how exactly the package (SummarizedExperiment) feeds data as an `assay` into a SummarizedExperiments object.
> 1 votes

---
Tags: rna-seq, deseq2
---


title: Repeating indexes of atoms in PDB file

text:
# Question
Title: Repeating indexes of atoms in PDB file
I have a PDB file in which the indexes of atoms are repeated (or better say refreshed) every 99999.

I am using MDtraj software to analyze this pdf file, so these repeating indexes cause some problems.

Does everyone encounter this problem and how can I fix it?

# Answer
> 1 votes
Inconsistencies like randomly repeated atom indices, maxed out indices at 99999, integer overflowing indices to the negatives (yup, there's a program that has this bug), bad spacing, jumbled order etc. in a PDB are tolerated differently in different programs. Some will carp at missing atoms or non-standard atom names, which are not formatting errors.

In terms of atom indices, these are used in `CONECT` entries, but several programs bond by proximity or by known topology —`CONECT` entries are absent for canonical amino acids. Repeated residue indices and chain letter repeats cause issues with `LINK` and often with peptide chain sequence.

However, in a multimodel PDB the numbers in each model are expected to be reset at every MODEL so if the atom indices reset after the end of a model this is 100% expected.

---
Tags: pdb
---


# Question
Title: What is the difference between Normalized Expression in EdgeR vs DESeq2?
I am trying to access the normalized expression in both edgeR and DESeq2, yet the results are different. Does anyone know why?

**How to get normalized expression using edgeR:**

```

y <- DGEList(...)
y <- calcNormFactors(y, method = "TMM")
v <- voom(y, design, plot=FALSE)

# normalized expression
head(v$E)

```

**How to get normalized expression using DESeq2:**

```

dds <- DESeqDataSetFromMatrix(...)
dds <- estimateSizeFactors(dds)

# normalized expression
normalized_counts <- counts(dds, normalized = TRUE)
head(normalized_counts)

```

Anyone know why the results are different? Aren't they both doing a version of log(cpm) under the hood? Or am I not comparing the right results from DESeq2?

# Answer
> 2 votes
First, and most obviously, two independent methods will never fully agree.

Second, while the DESeq2-related code in fact gets the normalized counts based on the size factors the edgeR/voom code extract the voom-transformed counts, but that is not the normalized counts itself. I would need to go back to the limma docs and paper to check what exactly the `E` slot stores, but I think (please edit for a definite answer) it's normalized counts after the voom transformation. In contrast, the edgeR normalized counts in the sense of what woul be most comparable to the DESeq2 normalized counts would be obtained with the `edgeR::cpm()` function. Check the manual. At no point of either the limma or edgeR manual it says that `v$E` should be used for anything outside of the linear modelling. It clearly recommends to use `cpm()`. Limma and edgeR are packages that have hundreds of forum questions and blog entries, many functions do exist and the packages change over time, so indeed it can be confusing what to do, but therefore I strongly recommend to consult the current manual at Bioconductor as the sole reliable reference on how to get values from the objects and how to run your analysis.

By the way, don't be surprised that the DESeq2 and edgeR normalized counts will be on different scales. DESeq2 gets the counts by dividing the raw counts by the size factors while edgeR first calculates naive per-million counts and then corrects these with the TMM factors calculated by `calcNormFactors`, so the magnitude of the counts can be quite different, simply by how it is calculated.

---
Tags: rna-seq, deseq2, normalization, limma, edger
---


# Question
Title: Extract chromosome arms' coordinates for specific reference
I would like to get chromosome arms' coordinates for genome reference `hs37d5` (ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/reference/phase2\_reference\_assembly\_sequence/), but I could not find them either online, or the directory listed above.

Is there a way to extract such coordinates from one of the available files (e.g., from the fasta)?

# Answer
The chromosome sizes can be retrieved with (to use for sanity check):

```
> ./faSize -veryDetailed -tab hs37d5.fa.gz
1   249250621   23970000    225280621   225280621   0
2   243199373   4994855 238204518   238204518   0
3   198022430   3225295 194797135   194797135   0
4   191154276   3492600 187661676   187661676   0
...

```

(first column is chromosome ID, second column is chromosome size).

Afterwards, the centromere coordinates are retrieved through (in this case I used hg19 since I'm working with `hs37d5`):

```
> curl -s "http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/cytoBand.txt.gz" \
  | gunzip -c | grep acen | sed 's/chr//g' | sort -n > hg19.centromeres

```

The coordinates for each chromosome arm can be computed from the chromosome size and centromere coordinates.

N.B. As per UCSC documentation, cytoband's start is 0-based, while end is 1-based. As a sanity check you might want to compare the end of the last cytoband for each chromosome, with the chromosome size (they should be the same, since end is 1-based).
> 1 votes

---
Tags: reference-genome, chromosomes
---


# Question
Title: Phylogenetic tree rooting in shotgun metagenomics
But I have some weeks fighting with this issue about phylogenetic tree building to use in a phyloseq object in order to calculate beta-diversity metrics that takes into account tree distance branches metrics.

I have one tree for Archea and another tree for bacteria (both unrooted), created from IQTREE from the protein alignments obtained by GTDB-Tk software. I tried to use the software Root-Digger, but only works with nucleotide alignments and not protein alignments. Do you know another software capable (working with protein alignments) to infer more objectively the root rather than define arbitrarily a root through ape R package?

---

Root-digger is described as "a non-reversible Markov model", which would be assumed to be a directional mutation matrix, thus T-\>A would have a different probability to A-\>T using Bayesian approach to branch-length estimation

IQ-TREE is an ultra-rapid method for non-parametric bootstrapping which can be essential for extremely large phylogenies.

# Answer
The answer really is that a root is determined using external evidence, rather than from the internal branches.

An algorithm such as root digger can be used to estimate the root placement based on a branch that for example might minimise the anomalous branch lengths in all other taxa. Excessive comparative branch-lengths could indicate an error in rooting the tree and this can be the basis of the estimate. However, it is an estimate.

R libraries such as `ape` are more likely to simply mid-point root the tree, which isn't really cool, but its easy to calculate.

The three approaches for rooting are either:

1. Using a known outgroup to the taxa of interest, which unequivocally predates their divergence.
2. External information such the evolution of key morphological features of the organism.
3. Published studies.

*Point 1* This can be used singly to find the root - a "root-finding tree" - and the tree of interest presented by constructing a second tree without this root in the alignment. For example the Archaea predate bacteria, one tree determines where the Archaea attach to taxa of interest. If the outgroup is very distinct this could cause artefacts.

*Point 2* Finding an outgroup much close to the taxa of interest using external information, for example gram-negative versus gram-positive bacteria, e.g. if all the bacteria in your analysis were one or the other and again do the two tree strategy (point 1).

*Point 3* Using large kingdom wide phylogenies from published studies to determine where the most ancestral branch in your taxa.

This is not an exhaustive list about how to use external information to determine the location of the root, but it is very widely recognised.

---

**Directional versus reversible matrices**

Directional matrix (e.g. most options in Beast) need to know where the root is, or they will 'guess it': ideally and this information can be feed into the algorithm when it is started. This is important if you are doing molecular dating.

In fact, you can do this even with reversible matrices (e.g. RAxML) but it makes no difference to the answer.

---

**Manual re-rooting**

If you have identified the root and the phylogeny you have does not represent that root correctly you can "re-root" the tree using the outgroup identified from point 1 to 3. A program such FigTree are cool for re-rooting. Re-rooting is perfectly cool providing you present the methodology that identifies the root taxa or taxon (singular).
> 2 votes

---
Tags: phylogenetics, phylogeny, metagenome, scaffold
---


# Question
Title: "Error in col2rgb(colors, alpha = TRUE) : invalid RGB specification" when running HeatmapAnnotation() (from ComplexHeatmap)
*This question was also asked on github*

I'm trying to run the function HeatmapAnnotation (a function in the ComplexHeatmap package). When I run it, I get the error "Error in col2rgb(colors, alpha = TRUE) : invalid RGB specification"

col2rgb() is not among the libraries I've loaded and it doesn't seem to be a method contained in HeatmapAnnotation.

Here's the code block for the section:

```
df <- data.frame(labid= combmetadata$Sample_Name,
             Diagnosis= combmetadata$Diagnosis,
             Tissue = combmetadata$Tissue,
             Tumor.Normal = combmetadata$T.N,
             study = combmetadata$batch)

 labid <- colnames(bValsselect)
 df <- left_join(as.data.frame(labid),df,by="labid")
 df <- df[,c(2,4,5)]
 df$study[df$study == "1" | df$study == "2"] <- "cgid"

 df$Diagnosis <- factor(df$Diagnosis, levels = c("Uterine endometrial","Ovarian endometrioid",
                                            "Ovarian mucinous",
                                            "Colorectal mucinous",
                                            "Pancreas mucinous",
                                            "Stomach mucinous"))

 ha_column = HeatmapAnnotation(df = df,
                          col = list(Diagnosis = c("Uterine endometrial" = "#DDCC77F",
                                          "Ovarian endometrioid" = "#0F7554",
                                          "Ovarian mucinous" = "#44AA99",
                                          "Colorectal mucinous" = "#882255",
                                          "Pancreas mucinous" = "#AA4499",
                                          "Stomach mucinous" = "#D695D0"),
                            Tumor.Normal = c("N" = "grey",
                                             "T" = "black"),
                            study = c("cgid" = "#332288",
                                      "tcga" = "#BD9AF7")))

hegt <- Heatmap(bValsselect,
            column_split=4,
            name = "Beta Value",
            col = colorRamp2(c(0,0.25,0.5),c("#00FFCC","#FFFFFF","#0099FF")),
            show_row_names = FALSE,
            top_annotation = ha_column,
            show_row_dend = F)

```

The code runs fine up until it reaches the end of the HeatmapAnnotation() call

I've tried to load col2rgb() and override the alpha setting. Can anyone give me a clue about what's happening and any suggestions on how to correct the problem? I'm running R 4.2. Thanks

# Answer
Looking at the documentation, it seems that `HeatmapAnnotation` expects a vector, whereas you are giving it a list. It is not clear on what acceptable color types are. It might even work with a list, but I'd suggest playing around with what you are supplying as your `col` parameter.

Note that `HeatmapAnnotation` has to do color mapping of some kind, the code is a little dense, but it does look like it calls out to `SingleAnnotation`, which in turn runs `ColorRamp2`, which does use `col2rgb`.

It also looks like `col2rgb` is called in `ColorRamp2`, which you use in the `Heatmap` call. How sure are you that your script isn't making it that far? It would help to see error output, possibly.

That said, when I take your isolated `ColorRamp2` and run it, it goes ok. Possibly pull out the `ColorRamp2` call and run it yourself in your environment to make sure nothing weird is happening.

Try moving the assignment of the colors outside the `Heatmap` call, occasionally supplying function calls as parameters like that causes issues.
> 1 votes

# Answer
The alphanumeric value for the first color assignment contains 7 alphanumeric characters instead of 6. Correcting this typo fixed the problem.
> 0 votes

---
Tags: r, bioconductor, heatmap, complexheatmap
---


# Question
Title: How to configure Cromwell to work with slurm?
I have a cluster that is set up with slurm and I seek to use Cromwell to dispatch jobs to the slurm scheduler. I followed the description explained on Cromwells Documentation Website. And created a configuration file called `cromwellconfig.conf`, however when I run it using the command `java -Dconfig.file=/lccs/home/logio9/Configure_Cromwell/cromwellconfig.conf \ -jar cromwell-77.jar run test-workflow.wdl` I get an error parsing the config file. I am confused as to where the error is, bellow is the contents of the configuration file:

```
include required(classpath("application"))

backend {
  default = SLURM
  providers {
    SLURM {
      actor-factory = "cromwell.backend.impl.sfs.config.ConfigBackendLifecycleActorFactory"
      config {
        runtime-attributes = """
        Int runtime_minutes = 600
        Int cpus = 2
        Int requested_memory_mb_per_core = 8000
        String queue = "short"
        String jobname = "testwork"
        String cwd = "/lccs/home/logio9/Configure_Cromwell"
        String err = "errors.output"
        String out = "ouput.out"
        """

        submit = """
            sbatch -J ${jobname} -D ${cwd} -o ${out} -e ${err} -t ${runtime_minutes} -p ${queue} \
            ${"-c " + cpus} \
        --mem-per-cpu ${requested_memory_mb_per_core} \
            --wrap "/bin/bash ${script}"
    """
    kill = "scancel ${job_id}"
        check-alive = "squeue -j ${job_id}"
        job-id-regex = "Submitted batch job (\\d+).*"
      }
    }
  }
}

```

This is the test wdl script I made, it copies a very large directory into another directory.

```
version 1.0

workflow Greetings {
    call SayHello {
    }
}

task SayHello {
    command <<<
        cp -R /lccs/home/logio9/Configure_Cromwell/tf1 /lccs/home/logio9/Configure_Cromwell/tf2

    >>>
    output {
        Array[String] out = read_lines(stdout())
    }
}

```

If anyone has had experience with this, please do share a bit on how exactly this configuration works.

# Answer
> 1 votes
Cromwell already defines the `cwd`, `err` and `out` runtime attributes by itself (it also defines a `job_name` and `script`). So when you define these variables in your config it causes a name conflict, since they already exist.

Cromwell will make a directory inside of the "cromwell-executions" directory (created when you run cromwell) for each job it runs. It will write the script, various stdout and stderr files, and the rc file to this directory. The intend is also for the job to run in that directory. Hence cromwell provides these runtime attributes to point to the appropriate locations within this directory so they can be used in your custom submission command.

You can try just removing these variables in the runtime-attributes section, while keeping them in the submit section.

---
Tags: wdl
---


# Question
Title: Handling failed blocks within scattered WDL workflow
I'm trying to adapt an existing workflow to make it more robust to failures. Below is a simplified example of the workflow.

```
workflow main {
   File filelist="list_of_input_files.txt"
   Array[String] files=read_lines(filelist)
   scatter(file in files) {
      call task1 { input: in=file }
      if(defined(task1.out)) { call task2 { input: in=task1.out } }
      if(defined(task2.out)) { call task2 { input: in=task2.out }
   }
}

task task1 {
   File in
   String out1 = basename(in,'.txt') + '_1.txt'
   command { if(...) { ./script1.sh ${in} > ${out1} } }
   output { File? out = out1 }
   runtime { ... }
}

task task2 {
   File in
   String out2 = basename(in,'.txt') + '_2.txt'
   command { if(...) {./script2.sh ${in} > ${out2} } }
   output { File? out = out2 }
   runtime { ... }
}

task task3 {
   File in
   String out3 = basename(in,'.txt') + '_3.txt'
   command { if(...) {./script3.sh ${in} > ${out3} } }
   output { File? out = out3 }
   runtime {...}
}

```

with my list\_of\_input\_files.txt looking like

```
firstfile.txt
secondfile.txt
thirdfile.txt
...

```

Obviously this is oversimplified, but I am trying to catch failed blocks by making the output of each task optional (using `File?`) and then having an `if(defined(taskX.out))` condition to call the next task. My reasoning is that if something fails (represented by the `if(...)` within the command block of each task), then the `out` file will never be created in that task and so won't exist as input for the next task and thus `if(defined(taskX.out))` will resolve as false and skip the next task (and therefore all the subsequent tasks).

With current testing, the first file in the list satisfies the `if(...)` conditions of all the tasks and succeeds in giving me the final output file `firstfile_1_2_3.txt`. However, I am getting a `Required file output 'secondfile_1.txt' does not exist'` error message because it fails the `if(...)` condition in `task1` and so the output file is never created.

From what I've read in the WDL v1.1 spec, optional `File` outputs are allowed and are undefined if the files they point to don't exist. However, after much confused googling, I've realised it's probably because I'm using Cromwell, which is stuck with WDL v1.0.

My question is whether this should be working in WDL v1.0 and, if not, is my logic correct and should this (psuedo-)code work with WDL v1.1? Is there a way to achieve what I'm trying to do using pure WDL v1.0?

# Answer
Optional outputs should work in WDL 1.0. Here is an example that runs without issue through cromwell (version 74).

Note that for the second call in the scatter the input file needs to be converted to from an optional to an non optional value. Since this only gets run if this value is defined, this conversion can safely be done using the `select_first` function (which returns the first defined value in a given list of optionals).

```
version 1.0

workflow test {
    input {
        File filelist = "files.txt"
    }
    Array[String] files = read_lines(filelist)

    scatter (file in  files) {
         call blah as b1 {input: in=file}
         if (defined(b1.out)) {
             call blah as b2 {input: in=select_first([b1.out])}
         }
    }
}

task blah {
    input {
        File in
    }
    String outPath = basename(in,'.txt') + '_blah.txt'

    command {
        if $(grep -q "yay" ~{in})
          then
             echo yay > ~{outPath}
        fi
    }

    output {
        File? out = outPath
    }
}

```

files.txt looks like this:

```
test1.txt
test2.txt

```

With test1.txt containing `yay` and test2.txt containing `nay`.
> 1 votes

---
Tags: wdl, workflow-management
---


# Question
Title: How to build advanced boxplots from datasets which columns contain only numbers?
I have got a fold-change dataset and I simply would like to plot something like this:

The problem is that I have only numbers in my columns and when I use a boxplot() function I can plot only boxes but not dots, jitter, etc.

I use this code to build a df and then proceed to a boxplot

```
A <- c(-1, 25, 3, 4, 50)
B <- c(23, -41, 32, 68, 26)
C <- c(63, 71, 22, -8, 2)
D <- c(16, -49, 2, 56, 6)
df <- data.frame(A, B, C, D)
print (df)
boxplot(df)

```

How can I plot something that looks like the above picture? Thanks

# Answer
> 2 votes
To coerce your data into a form interpretable by `ggplot2`, you can use the `gather` function from `tidyr`.

`gather` has been superseded by `pivot_longer`.

```
library(tidyverse)
df <- tibble(A = c(-1, 25, 3, 4, 50),
             B = c(23, -41, 32, 68, 26),
             C = c(63, 71, 22, -8, 2),
             D = c(16, -49, 2, 56, 6))

tidyr::pivot_longer(data = df, cols = c(A,B,C,D)) |>
  ggplot(mapping = aes(x = name, y = value)) +
  geom_boxplot() +
  geom_dotplot(binaxis = "y", stackdir = "center") +
  theme_classic()

```

---
Tags: r
---


# Question
Title: Transfer learning for a Convolutional NN for Recursion dataset
As a practice data analysis I am trying to train a convolutional neural network (CNN) on some cellular images made publicly available on Kaggle in 2019 by the company Recursion via supervised deep learning.

I want to increase the predictive power using transfer learning at the lower layers of my CNN preferably from a public imaging model.

Are there any appropriate models available?

# Answer
```
learn = cnn_learner(data,
                    models.resnet50,
                    metrics=error_rate)

```

resnet50 is one of many pre-trained 2D CNN computer vision models which use transfer learning. If you've classifying a common image such as a given mammal you're sorted because they'll be loads of pictures of that mammal all over resnet50. However, 'a cell' I dunno.

When the computer vision challenges were in full swing between the tech giants each put their trained models onto publication release and are easily accessible e.g. via keras, fastai etc... I'd have to look them all up and each has different architecture, they are serious architecture like 50 layer NN. The competition was eventually stopped because computer vision had a lower error rate than human vision, at which point its pointless judging between DL algorithms because we can't say what it is how can we judge an algorithm?

Transfer learning places CNN at significant advantage over RNN.

The other possibility is google images

```
pip install google_images_download

```

I don't think its on conda. Then see if you get lucky with the following

```
from google_images_download import google_images_download
response = google_images_download.googleimagesdownload()
arguments = {"keywords":"cells,organelles, membranes","limit":60,"print_urls":False}
paths = response.download(arguments)
print(paths)

```

Please note I don't do 2D CNN,i.e. images.
> 1 votes

---
Tags: machine-learning, deep-learning
---


# Question
Title: How will Seurat handle pre-normalized and pre-scaled data?
I want to use data from the following dataset: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE84371

The data has been TPM normalized, which is not ideal for clustering but I have to work with what I have. TPM is a bad normalization method and it should not be used for these analyses because its laden with a lot of assumptions.

Presumably it has already been scaled, ***can I simply skip these steps and then run FindVariableFeatures() followed by the rest of the pipeline from the clustering tutorial or do the methods used assume Seurat specific preprocessing and using TPM will screw them up?***

Unfortunately in this case TPM is not reversable. https://btep.ccr.cancer.gov/question/faq/what-is-the-difference-between-rpkm-fpkm-and-tpm/#:~:text=Here's%20how%20you%20calculate%20TPM,divide%20this%20number%20by%201%2C000%2C000.

TPM is calculated as a function of the kilobases of the gene in the reference genome they used and the reads per gene across all samples (aka all cells). They have not specified which reference genome was used in the main article, and they do not supply the read counts in the metadata or the main file. So in this case the TPM is irreversible.

I think I can get raw reads from SRA and do alignment myself. It's annoying that raw reads were not provided but as someone who's lab entirely does computational analysis on open datasets I can tell you that this is not the most egregious thing I've ever seen. As I said before it would be easier for me to get the fastq files and do my own alignment which is still a ton of work the main authors could have prevented by simply posting the read counts.

Clustering Tutorial: https://satijalab.org/seurat/articles/pbmc3k\_tutorial.html

---

TPM means transcripts per million. It is obtained by dividing the read counts by the length of each gene in kilobases, resulting in reads per kilobase (RPK). All the RPK values in a sample are counted and then divided by 1,000,000.

# Answer
I don't do transcriptome analysis, it ain't my thing, however I do understand statistical analysis as well as the underlying issue regarding the public availability of molecular data ... I agree with the OP its not ideal. However, yes the OP can continue with 'clustering', personally I definitely prefer it if the transformation was removed, or analysed from the raw data.

**Summary** The authors of the data set are obliged to provide the raw data, presuming its published. Why? This is because whilst the OP can perform the analysis - it is very likely the OP will be asked to repeat the analysis on the untransformed data (again note this ain't my area, so 'very likely' might be 'possibly' in this field). At this point the original authors have to declare their untransformed data, if this was a formal request to the OP's analysis.

Please read the caveats and discussion.

1. Firstly could the OP reverse the normalisation? Its just a transformation applied to all the data which in their analysis resulted in a normal distribution, I don't know TPM but if this was log-transform you would supply 'un-log' it, square-root then square it, reciprocal do another reciprocal will 'de-transformation'. Can TPM be undone? I agree it ain't ideal, from several perspectives but it gets things moving quickly.
2. Ask the authors for the raw data - this is what I would do. They are obliged to provide it if the data is published. Virtually all journals have this policy in place, however I agree that ASM remain an outlier in this respect (who are significant publisher). Any other journal will return to the authors and point out their 'check-boxes' when they submitted the manuscript. It is etiquette to go to the authors first and provide sufficient chance for them to respond, prior step 3 (below) - giving them the absolute maximum benefit of the doubt.
3. What you do if the authors don't respond? The OP would then contact the journal explaining: a) why you need the data, i.e. the analysis performed; b) explain the analytical significance of the raw data, i.e. arguments the OP has raise 'not ideal' and general concerns raised here. I do have experience in this regard and the authors will oblige if asked by the journal.
4. Searaut's cluster analysis certainly used tSNE ... I believe it uses dimensionality reduction such principle component analysis (PCA) - I don't know which they use - and then tSNE analysis (which I think you are calling 'clustering\*'). This is a very robust 'clustering'\* (please see footnote), there is no question this is as good as it gets. If that is singly your goal I agree that dimensionality reduction calculations don't need transformations, transformations are needed for parametric statistics and other areas of model fitting.

**The problems with transformation**

I agree the authors should not have done this because:

1. Transformations are not perfect, if its wrong - you're stuck.
2. Not only do you not need to transform the data, its better you don't - if untransformed data results in nice 'clusters' thats the preferred result. Dimensionality reduction + tSNE gives great results thats well known.

**Corroborating the OP's concern** Dimensionality reduction calculations do not need a normal distribution because its sectioning up the observable variance to maximise the classification of the data. My opinion - in statical analysis - is if you need to transform, THEN perform dimensionality reduction (PCA-stuff) THEN tSNE, that data has been heavily manipulated. It doesn't mean to say the conclusions are not valid but its been put through a lot and would STRONGLY SUGGEST THE RAW DATA HAS NOT GIVEN GOOD RESULTS. Thus, personally I would be cautious, why should be OP be put under unnecessary caution?

It is important to note this is not my field and the field will have its own standards. If transformation is routine for all statical calculations (except non-parametric where it may not be valid) then authors of the data will argue they are correct and no further data needs to be released. Personally, I really doubt it, but we'll see what others say.

**Footnote** * clustering is a statistical method usually resulting in a tree, it is generally avoided because there are so many different clustering methods giving different results and no agreement on which is the 'optimal' method. I think when the OP say 'clustering' they refer to clusters following a tSNE analysis, which will be a bivariate plot.

---

**From the comments** Failure to declare a reference genome is a formal omission that has to be rectified and the journal is obliged to do this. Thus the OP would simply contact the journal, they may tell the OP to contact the authors and when the authors respond the response is passed back to the journal. If its an online journal you can simply leave a comment in the 'comments' section: the authors will answer. In this instance officiating via the journal is preferred to ensure the information is verified by all parties - I presume the OP does not know the authors. From the comments I get the gist the OP is working via supervision, and would therefore just notify their supervisor who should officiate on OP's behalf (they might then delegate it back to the OP - but at least the supervisor is formally aware of this).

In context, stuff happens and there is formal process to ensure parity is achieved its simply patience in working through it.

---

**Down vote** In response to the downvote ... unpopularity happens in data set transparency its part of the process ... I was (at least in part) attributed to changing journal group policy on data availability, stated from within the journal. Maybe they joking - but they certainly all knew about what was going on. The specific issue was a prolonged dialogue over access to a specific data set. It does demonstrate that most journals are concerned if an author does not cooperate with data transparency, in fact the journal in question explicitly said this and they were not joking.
> 1 votes

# Answer
I notice that the paper points to the Broad website for data accession; have you looked there?

https://singlecell.broadinstitute.org/single\_cell/study/SCP1/-single-nucleus-rna-seq-of-cell-diversity-in-the-adult-mouse-hippocampus-snuc-seq

For specifics on the normalisation method, you could look at their associated github repository:

https://github.com/yinqingl/nucseq\_analysis

Given that the paper was published in 2016, they are TPM counts, and you have access to the raw reads, my recommendation would be to regenerate gene count tables from those raw reads. I expect it will be more time consuming, less accurate, and require more effort to track down the gene count table than to create one from scratch. You can use Salmon Alevin to do this fairly quickly:

https://salmon.readthedocs.io/en/latest/alevin.html

\[bearing in mind that there may be some footwork required to work out how to get Salmon to interpret the barcoding method that was used\]
> 0 votes

---
Tags: rna-seq, scrnaseq, seurat
---


# Question
Title: What does this accession NCBI code mean: 6MWN_B?
According to this article, accession codes should consist from a combination of uppercase letters following a combination of digits. If this is a *RefSeq*, it can have a prefix as a combination of uppercase letters with underscore.

But this accession code **6MWN\_B** does not consist with all aforementioned. What is it? I could not find an answer via Google.

# Answer
> 3 votes
In general, an *accession* is a value that points to (provides *access* to) a single database entry. Different databases (GenBank, RefSeq, Ensembl, SwissProt, and so on) have different conventions for accessions.

The accession in question points to an RNA sequence associated with a Protein Data Bank (PDB) record: "Crystal structure of hepatitis A virus IRES domain V in complex with Fab HAVx". In the PDB record, you can click the "sequence" tab and select chain "B" from the drop-down, and you will see an RNA sequence (ACGU) matching the DNA sequence (ACGT) in the NCBI record to which you linked.

I'm not sure what the rules or conventions for PDB accessions are, but they are clearly different from RefSeq accessions.

---
Tags: phylogenetics, genbank
---


title: ATAC peak anatations

text:
# Question
Title: ATAC peak anatations
I am new in ATAC analysis.

I have a problem about how to transfer ATAC peak to gene symbol.

Here is the ATAC-seq: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM3271041

The `GSM3271041_ATAC_sciCAR_A549_peak.txt.gz` file contains five columns : `(id, peak, chr, start, end)`. What I am doing is to use `bedtools` to transfer locations in chromosome to get gene symbols.

Here is an example from `GSM3271041_ATAC_sciCAR_A549_peak.txt.gz`

```
id,peak,chr,start,end
9,1-68591-69091,1,68591,69091
10,1-91105-91605,1,91105,91605

```

I use bedtools to get the gene symbol of `id=9` , that is `OR4F5`.

I want to use these ''symbol-transfered'' ATAC-seq to integrate with RNA-seq which has gene symbols alreadly.

Am I right to do so ?

# Answer
ATAC-seq data is fundamentally different from RNA-seq data because they measure different things. RNA-seq data is typically mapped from raw sequencing reads to annotated gene transcripts. This has a nice 1:1 correspondence because RNA is transcribed from genes.

ATAC-seq data, on the other hand, is used to measure accessible regions of the genome. But accessible regions don't come from genes the same way that RNA comes from genes. Accessible regions are scattered throughout the genome and most lie in intergenic regions where there is no annotation for a gene. Because accessible regions don't have the same biological connection to genes that RNA transcripts do, you need to be more precise in how you want to pair up accessible regions and genes.

Are you looking for *cis*-regulatory elements (CREs) that regulate nearby genes? Those will often be found within accessible regions, so you'd need additional data that connects an individual regulatory element to that gene.

Are you looking to create a catalogue of accessible regions within a certain cell type? Then you might want to look at ENCODE's candidate CREs database to see how your ATAC-seq peaks compares. From there, you'll be able to find additional information about those peaks, such as whether other functional genomics data identifies these regions as active regulatory elements, which gene(s) these elements regulate, and more.

There are different options depending on what exactly you're trying to do, but there is no one way to transfer ATAC-seq peaks to gene IDs because there's no biological question there.
> 3 votes

---
Tags: atac-seq
---


# Question
Title: How many single copy gene are needed to construct the phylogenetic tree?
I ran the orthomcl program with genomes of the whole streptomyces genus(from ncbi) as input, but I only got 22 orthologous groups(streptomyces's genome size vary from 2M to 15M, which shows the great difference among them). And the analyst from the bioinformatics company told me that my result wasn't reasonable if I try to construct the phylogenetic tree with the 22 orthologous groups I got. He also told me, I'd better use almost 1000 orthologous groups to build the tree. My question is, how to know my orthologous groups are enough to construct a tree?

The program was run during my summer vacation(a long time ago...), but I still can't realize why "1000 orthologous groups are needed". And the analyst form company just tell me that it was the experience number. But when I read some article about comparative genomics, I see some papers are used just about 10 single-copy genes to build the tree, why???

# Answer
The one thing that strikes me is that OrthoMCL was written under David Roos supervision, which really means it was written with Toxoplasma gondii in mind, i.e. a eukaryotic parasite with notable levels of gene duplication. They tested it against Plasmodium falciparum, but that isn't surprising because both are members of the Apicomplexa. Personally I would suggest looking at a tool more orientated towards bacteria, albeit I'd like to see the output of OrthoMCL.

Prokka is a tool that might give better results, here, which is designed for bacteria, but personally I would extract the orthologes/homologues directly if the primary purpose is to contruct a phylogeny, i.e. I would parse annotated genes which are not "hypothetical".
> 1 votes

# Answer
I would say, 120 markers is enough for phylogenetic analysis. see https://gtdb.ecogenomic.org/
> 1 votes

---
Tags: phylogenetics, phylogeny, bacteria
---


# Question
Title: Where to download a table with ICD-9-CM codes?
I am looking for a simple table with all ICD9 codes (For example: Diagnosis Code 414.0 -\> Coronary atherosclerosis). I can only find weird PDF or other kinds of text documents. Or web forms, where I have to post a certain code to get the description.

## Here is what I have found so far.

http://www.icd9data.com/2015/Volume1/default.html

https://www.aapc.com/codes/icd9-codes-vol3-range

https://www.cms.gov/Medicare/Coding/ICD9ProviderDiagnosticCodes/codes

https://www2.gov.bc.ca/gov/content/health/practitioner-professional-resources/msp/physicians/diagnostic-code-descriptions-icd-9

https://ftp.cdc.gov/pub/Health\_Statistics/NCHS/Publications/

# Answer
The zip files found within https://www.cms.gov/Medicare/Coding/ICD9ProviderDiagnosticCodes/codes contain text files and Excel spreadsheets that you can open with Excel, R or Python depending on how you want to use this data.
> 1 votes

---
Tags: icd-codes
---


# Question
Title: Tree Building Algorithm that treats gaps as deletions
I'm part of a nanopore sequencing experiment that will sequence several generations of viruses. The intent is to perform directed evolution by putting selective pressure on these viruses and tracking the various mutations that occur. See here for the paper this methodology is based on.

I've been working on a software package that can compile consensus sequences from longread nanopore sequencing (see here). However, there's a wrinkle with generating phylogenetic trees from this type of data. It appears that most tree building software ignores indels. This makes sense, as gaps could easily indicate a lack of coverage from illumina sequencing methods or similar issues rather than betray the existence of an insertion/deletion. However, the nature of our sequencing method is such that a gaps really should be evaluated as mutational changes that we want to keep track of in a phylogenetic tree.

We've explored using NextStrain to generate our trees and visuals, but the tree building algorithms they employ too suffer from this gap assumption (see the issue I raised here). These include IQTree, RAxML and PhyML. I have not had much luck finding an algorithm or software that treats gaps as indels when building a phylogenetic tree. Does anyone know of such an algorithm or package I could employ? I figure I could take a stab at writing my own algorithm, but consider that a Plan B resort. No sense re-inventing the wheel and all that.

# Answer
I need to do a more thorough testing to determine just how effective this answer may be, but I think I've put together a workflow that is close enough to an appropriate answer that I'd post it here.

In my research I came across gappy, (GitHub here), a piece of software that identifies potentially relevant gaps in a phylogenetic analysis and returns a PHYLIP-formatted file that can be run through most tree-building programs. Gappy accounts for indels, but not substitutions, so the output would need to be combined with another analysis to generate the tree I want.

Fortunately, IQTree allows for a partitioned analysis of mixed data (see here), so by feeding the original alignment AND the gappy output through IQTree you should be able to generate a tree that accounts for both substitutions and indels.

Again, more thorough testing is required, but I tested this possibility using an MSA example taken from this website. Saving the FASTA-formatted example as `alignment.fasta`, and running that file through gappy2, you get an output file called `GAPPY2_alignment.fasta_2-inf.phy`. Thus if you create a NEXUS-formatted file like the following:

```
#nexus
begin sets;
    charset part1 = alignment.fasta: *;
    charset part2 = GAPPY2_alignment.fasta_2-inf.phy: *;
end;

```

and run it through IQTree using the command `iqtree -p example.nex`, you get the following tree:

```

+---ENA|CAA23748|CAA23748.1
|
+-------------ENA|CAA24095|CAA24095.1
|
|              +-------------------------------------------ENA|BAA20512|BAA20512.1
+--------------|
               +-----------------------ENA|CAA28435|CAA28435.1

```

Now, disclaimer: this is the exact same tree you get when you run `alignment.fasta` through IQTree, hence why more thorough testing is required. The sample size for this example is quite small, and it's therefore not surprising that accounting for indels does not change the outcome in this scenario. However, I think this serves as a solid baseline workflow, where future updates will be more tweaks than complete overhauls, if that makes sense.

Let me know if you see any glaring flaws in my logic, but I think I'm going to close this question. I may provide edits to this answer in the future if significant tweaks/overhauls wind up being done.
> 1 votes

---
Tags: phylogenetics, nanopore, long-reads, indel
---


title: Violin plots appear as vertical lines

text:
# Question
Title: Violin plots appear as vertical lines
I want to perform a violin plot from a Seurat object and I rather have vertical lines. How can I stop this and rather have a standard Violin plot?

code

```
Seurat::VlnPlot(object, features=c("LTA", "LTB"), group.by = "group",
    pt.size = 0)

```

# Answer
> 3 votes
I believe the reason is simply LTA and LTB not being expressed (or their expression not being detected due to technical limitations). If you would set the `pt.size` argument to something else than 0, you would see a big blob of dots at point 0.0 (depending on the number of cells/data points).

---
Tags: seurat, ggplot2
---


# Question
Title: Normalization methods to combine scRNA-seq experiments with different sequencing depths
I am training a classifier to identify a cell type in a particular state of activity using scRNA-seq. There is a large variation in the sequencing depth (reads average per cell) of the testing data (data I will be classifying not training or validating the classifier with) and my supervisor suggested binarizing the genes. I was wondering if there are other options I should consider. I was wondering if there are particular normalization methods for ML classifiers which would be best for this task.

# Answer
> 1 votes
Coincidentally, a recent preprint out of Lior Pacther's lab covers some of these details around scRNA-seq depth normalization.

They find that a "PFlog1pPF" normalization method is able to stabilize variance across genes while still reducing the impact of total cell depth on that variance. If I've read this preprint correctly, "PFlog1pPF" is essentially a transformation of the form $\frac{1}{N\_c} \log \left(1 + \frac{k\_g}{N\_c} \right)$ where $N\_c$ is the number of reads for cell $c$ and $k\_g$ is the number of reads for feature $g$.

This preprint is mostly dedicated towards differential expression and identifying marker genes in certain cell populations, but that may still help you in your application.

Using a transformation like this (or others that are discussed in this preprint and its references) gives you more granularity and sensitivity for highly expressed genes than what binarizing the data will give you. And it has additional benefits like stabilized variance, which can cause problems in some downstream applications if you don't keep it under control. But starting with a binarized input like your supervisor suggested might be a simple first attempt to see if your idea has some potential in the first place.

# Answer
> 0 votes
I don't think you need to complicate the idea of normalisation by introducing machine learning classifiers as a necessary component. Normalisation is common when comparing different datasets for all differential analysis.

Given that you have single cell data, have a look at integration techniques in the Seurat workflows. These include normalisation as a component of integration:

https://satijalab.org/seurat/articles/integration\_rpca.html

---
Tags: scrnaseq, transcriptome, normalization
---


# Question
Title: What is the best way to perform batch ajustment between two experiments?
I have duplicate drug efficacy experiments performed six months apart. I would like to know what is the right way to perform statistical analysis to understand fold changes after treatment with a drug.

Say I have

```
ctrl1       tube1       ctrl2       tube2
7,732637722 10,71639307 5,454318175 8,052899428
6,69630905  8,762916793 7,441582994 9,855106201
10,35624396 11,73716043 8,976524001 11,0975511
4,29150392  4,848632578 5,901132154 7,918630661

```

and the same experiment but 6 months after (numbers are slightly different)

Could I take the mean between samples divided by two can work out? If not - why?

I would like to have statistically significant values for futher analysis. Thanks.

# Answer
> 0 votes
The answer is simply you use both sets of values in a one-way paired ANOVA. Note that certain criteria need to be in place for ANOVA, notably the Fmax test.

Looking at the data 2-way ANOVA might be needed and that is more complicated, or else 2 one-way ANOVAs: I don't understand 'ctl2' and 'tube2' data. Is this a single or multiple drugs?

Wilcoxon signed-rank test is the first choice test because there are no prior assumptions required, which replaces both the T-test and the 1-way ANOVA.

Whatever you do, taking the mean and performing a T-test is not a good idea. You asked why ... because the variation between the values is critical to the significance of the test and is simply demonstrated:

Tube 1 value

```
tube1 time X = 11
tube1 time y = 9
tube 1 mean = 10

```

Tube 2 value

```
tube 2 time x = 19
tube 2 time y = 1
tube 2 mean = 10

```

Thus both would appear to have no change - they are both 10 - and yet the variation between them means they observe very different behaviours. Thats why.

---
Tags: statistics, batch-effects
---


# Question
Title: how to create "sample file" for the qAlign() function after trimming the reads in R
I'm an absolute beginner trying to solve this question "Align the trimmed and untrimmed reads using QuasR and plot alignment statistics, did the trimming improve alignments?"

I did trim the reads using the preprocessReads() function, now I want to align the output but the qAlign() only takes a samples file and genome file. however, my trimmed reads are two separate files how can I align them.

I asked google but I didn't find the answer

# Answer
The QuasR vignette tells you how to run it. The samples file should contain a list of your samples and fastq files.
> 1 votes

---
Tags: r, rna-seq, sequence-alignment, homework
---


# Question
Title: Association analysis for methylation beta value and phenotypes in R
I am trying to learn how to do association analysis of methylation array data with numerical trait like body weight, eGFR. As I have seen some of the publications, methylation data association with eGFR has been showed through regression analysis but I am not sure how to run regression for so many CpG probes in R.

Can you have look on this below publication and could you help with R code for regression analysis?

https://www.nature.com/articles/s41467-017-01297-7.pdf

I do have methylation beta value from RnBeads but I am not sure how to read methylation data (rows/column-CpG/SampleID) and then how to read phenotype data file and finally how to make single sheet to run regression analysis?

# Answer
> -1 votes
*Answer from @user4939328, converted from comment. Please improve this answer if you can:*

Have a look at this method:

https://academic.oup.com/bioinformatics/article/34/15/2657/4939328

*\[Wondering why this was converted from a comment? See the StackExchange rules about comments here\]*

---
Tags: r, bioconductor, statistics, epigenetics
---


# Question
Title: VCF or BAM file for raw data of gene test?
My friend has a VUS (Polymicrogyria with or without VEDS) mutation that was found in her whole-exome sequencing with respect to the phenotype given that time. At present, her doctors have more information on the complexity and are looking for more explainable mutation present in herself, if any. We asked the company and they're willing to give the raw data.

**Which file format should we opt for?** They've given the option of BAM & VCF.

# Answer
First of all - I would advise you against concluding anything out of your own analysis of the data, you need to consult someone with appropriate training!

Can you get both? The `.bam` file is a format that stores information about sequencing reads (usually mapped, but not always). You can even extract the original sequencing information from bam files. This is practically the rawest information you can have. The `.vcf` file is a variant file that gives you some information about inferred variants, but there are plenty of bits of information that will be missing. The basic information regarding sequencing depth and quality of variant calling is still included and therefore most likely sufficient for your physician to see any other variants potentially relevant to your friend's condition.

I will repeat once again - don't conclude anything on your own. It takes years to gain expertise in analysing genemics data and there is very high risk of misinterpretations if not analysed with care.
> 1 votes

---
Tags: sam, vcf, gene, mutations, genetics
---


# Question
Title: How can you put a bash array into a wdl variable?
Given two separate wdl arrays, my command section will do some modifications to each array and place the output into two separate bash arrays. Finally, I want to output the two separate bash arrays in the output section of the wdl task to be used in other tasks within the workflow. The below code is an example of what I would like to do:

```
task AddA {
    input {
        Array[String] wdllistX = [ "True", "False"]
        Array[String] wdllistY = [ "Good", "Food"]
    }

    command <<<
        ###We add the letter 'a' at the end of each string in each wdl list and finally get the
      bash arrays:

        bashX = ("Truea" Falsea")
        bashY = ("Gooda" "Fooda")
    >>>
    output {
        Array[String] finX = bashX
        Array[String] finY = bashY
    }
}

```

The above example produces an error. It seems that I might have to convert the bash arrays back into wdl arrays in order to declare them as outputs in the output section, however, I am not sure how to do this. May someone help me with this? Moreover, a general explanation of how you can capture the output of specific things from the command section into the output section would be helpful aswell.

# Answer
When a WDL task is run, the execution engine takes the command block and evaluates it as a template, rendering any wdl expressions it contains into the final script that is executed, **BEFORE** it is actually executed.

The actual execution environment is generally isolated from the execution engine or from other tasks. Information is passed between tasks by reading from or passing around file based outputs and not from accessing the bash environment directly, which the engine does not have access to.

So, in order to achieve what you are trying to do, you first need to get the values into a file format, and then you need to read those values using an engine function provided by `WDL`. Generally, all outputs need to be in a file if you want the engine to be able to interact with them. Please note that `stdout` and `stderr` are considered "files" and are accessible with `stdout()` and `stderr()` engine functions.

## Option 1: Construct json files and use read\_json

```
task AddA {
    input {
        Array[String] wdllistX = [ "True", "False"]
        Array[String] wdllistY = [ "Good", "Food"]
    }

    command <<<
        ###We add the letter 'a' at the end of each string in each wdl list and finally get the
      bash arrays:
        echo '["Truea","Falsea"]' > x.json
        echo '["Gooda","Fooda"]' > y.json
    >>>
    output {
        Array[String] finX = read_json("x.json")
        Array[String] finY = read_json("y.json")
    }
}

```

## Option 2: use read\_lines

```
task AddA {
    input {
        Array[String] wdllistX = [ "True", "False"]
        Array[String] wdllistY = [ "Good", "Food"]
    }

    command <<<
        ###We add the letter 'a' at the end of each string in each wdl list and finally get the
      bash arrays:

        echo "Truea" >> x.txt
        echo "Falsea" >> x.txt
        echo "Gooda" >> y.txt
        echo "Fooda" >> y.txt
    >>>
    output {
        Array[String] finX = read_lines("x.txt")
        Array[String] finY = read_lines("y.txt")
    }
}

```
> 4 votes

---
Tags: wdl
---


# Question
Title: Wrapping R scripts in WDL and run in docker
I am trying to run an R script as a task in WDL in a docker image. Since this is my first time working with these, I am stuck in something very simple (and in the initial stage) and am hoping somebody here can point me in the right direction. Here is my dockerfile.

```
FROM rocker/r-base:latest
RUN Rscript ...r-docker/install_packages.R

```

Here is my WDL script

```
version 1.0

workflow Rpipe {
   input {
     String name
   }
  call PrintNameR{
        input:
        name = name
}
  output {
     File outfile = PrintNameR.output_R
  }
}

task PrintNameR {
      input {
           String name
}

command<<<
    R --no-save --args name <<Rscript
    args <- commandArgs(trailingOnly = TRUE)
    for (i in args){
    print(i)}
    Rscript
  >>>

runtime {
            docker: "firsttrydocker"
       }
output {
     File output_R = stdout()
  }
}

```

I validated, generated an input json file - edited it to include the string I want "Smith", and ran it.

```
java -jar ../cromwell/womtool-77.jar validate rtry1.wdl
java -jar ..cromwell/womtool-77.jar inputs rtry1.wdl > rtry1.json
java -jar ..cromwell/cromwell-77.jar run rtry1.wdl --inputs rtry1.json

```

I get an output:

```
> args <- commandArgs(trailingOnly = TRUE)
> for (i in args){
+ print(i)}
[1] "name"

```

I get "name", I wanted to get "Smith" because that's what my input is.

What am I doing wrong? Clearly, it seems that I am not passing the WDL task input to the WDL task command -in the Rscript. What is the proper way to do it?

# Answer
In order to reference a WDL input within a command block you need to use an expression placeholder. Expression placeholders have the form `~{*wdl expression*}`, where the inner part can be an valid WDL expression (ie reference to an input declaration)

When an engine sees a `String` or a `command` block, it first looks to see if there are any of these expression placeholders within it. If there are, it will evaluate the inner expression, and replace the entire placeholder `~{...}` with the resulting value.

You can see how this can be applied below

```
String first_name = "John"
String last_name = "Smith"
Int age = 45

String full_name = "~{first_name} ~{last_name} # equals "John Smith"
String name_and_age = "John is ~{age}" # equals "John is 45

# Any valid WDL expression can go within an expression placeholder
String complex = "~{first_name + last_name + age}" # equals "JohnSmith45"

```

A key thing to note is that WDL command templating and expression evaluation is completely separate from the actual execution of the command. An execution engine will evaluate the command template, swap the expression placeholders with their evaluated values, and *THEN*, execute the script.

This is relevant for your use case since you are using heredoc to directly define the rscript you are runnning. If we briefly look at a python example, you can see what I mean:

```
task hello_world {
  input {
   String name = "john"
  }

  command <<<
   python <<CODE
     # this is valid, because once the placeholder evaluates
     # this command becomes name = "john", which is valid python syntax
     # We did not have to pass ~{name} as an arg to the python command
     name = "~{name}"
     print(name)
   CODE
  >>>
}

```

---

Now getting to your problem specifically. You are trying to reference the `name` input value without using any sort of expression placeholder. To the execution engine, the `name` value looks like any other string. To swap the `name` in the command for the value passed into the task, you simply need to wrap the variable name in an expression placeholder: `~{name}`.

task PrintNameR {
  input {
    String name
  }

  command<<<
    # I would wrap this in quotes to prevent any side effects of
    # shell expansion
    R --no-save --args '~{name}' <<Rscript
    args <- commandArgs(trailingOnly = TRUE)
    for (i in args){
    print(i)}
    Rscript
  >>>

  runtime {
            docker: "firsttrydocker"
  }

  output {
     File output_R = stdout()
  }
}
```
</code></pre>
> 2 votes

# Answer
I figured I had to pass ~{sep=" " name} as an argument

```
command<<<
R --no-save --args ~{sep=" " name} <<Rscript
args <- commandArgs(trailingOnly = TRUE)
for (i in args){
print(i)}
Rscript

```

> > >
> 4 votes

---
Tags: r, wdl, docker
---


# Question
Title: Python; Question about the append (modifying lists): why modifying a list with append changes others?
I have a list named dict\_eg1 whose values are lists.

```
dict_eg1 ={
    "A":[1,2,3],
    "B":[11,22,33],}
val1=dict_eg1["A"]
val1.append (100)#modify the list by using append.
print (dict_eg1)

```

{'A': \[1, 2, 3, 100\], 'B': \[11, 22, 33\]}

In this case, dict\_eg1 changed.

However,

```
dict_eg1 ={
    "A":[1,2,3],
    "B":[11,22,33],}
val1=dict_eg1["A"]
val1 = [1, 2, 3, 100]#modify the list without using append.
dict_eg1

```

{'A': \[1, 2, 3\], 'B': \[11, 22, 33\]} In this case, dict\_eg1 did not change.

Can anyone tell me why these two results were different??

# Answer
> 1 votes
Welcome to Bioinformatics.

First of all, it should not be asked in this community since it is purely a question about python programming. I suggest that next time, please head to StackOverflow as you will likely receive more help from there.

You also look new to programming. This is more of a logical question. Let's break it up to understand what is happening in your first example:

```
# you define your dictionary containing two keys along with their values.
dict_eg1 = {
             "A":[1,2,3],    # Key 'A' with a value as a list [1,2,3]
             "B":[11,22,33]  # Key 'B' with a value as a list [11,22,33]
}

# you now define a new variable 'val1' which is
# the value of key 'A' and that is the list [1,2,3]
val1=dict_eg1["A"]

# as you have mentioned, doing the below command
# will append the value to the list of key 'A'
# in other words, you can also write it like so
# dict_eg1["A"].append(100)
# that's why when you make a change to 'val1'
# the value of the key 'A' of 'dict_eg1' will change
val1.append (100)

```

As for the second example (I will skip the same steps):

```
# ...

# you store the value of key 'A' in 'val1'
val1=dict_eg1["A"]

# however, in the below command, you actually
# !redefine! the value of 'val1' to a new list
# that is totally independent of your 'dict_eg1'
# In this case, 'dict_eg1' is not affected by
# any changes made to 'val1'
val1 = [1, 2, 3, 100]

```

Just to sum it up, when you use `val1=dict_eg1["A"]`, imagine that you link the value of your key `A` to `val1` and any changes made to `val1` will affect `dict_eg1["A"]`. However, ***redefining*** `val1` will disconnect that link.

---
Tags: python
---


# Question
Title: ggtree：How to keep leading zero ‘0’ in the tiplabel？
I am using `ggtree` to draw a phylogenetic tree.

Some taxa IDs have a 0 at the beginning.

However, when I plot the tiplabel the leading zero is automatically dropped. How do I retain the leading zero?

Thanks.

The code and tree file and metadata file are attached. I would appreciate any help.

```
setwd('G:\\Experiments\\Paper\\Torula\\06_iqtree')
dir()

# Load required packages
library(xlsx)
library(ggtree)
library(treeio)
library(ggplot2)

# Read tree
tre <- read.iqtree("concatenate_iqtree2.treefile")
tre_tbl <- fortify(tre)
# Read associated information table
metadf <- read.xlsx('G:\\Experiments\\Paper\\Torula\\00_data\\Torula_taxa.xlsx', sheetIndex = 1)

# Plot tree
ggtree(tre) %<+% metadf +
    # our own taxa darkblue
    # type strain bold
    geom_tiplab(aes(subset = (New == 1 & Type == 1),
                    label = paste0('bolditalic(', Genus, ')~bolditalic(', Epithet, ')~bold(', Collection,')~bold(', Number,')')),
                color = 'darkblue', parse = T) +
    geom_tiplab(aes(subset = (New == 0 & Type == 1),
                    label = paste0('bolditalic(', Genus, ')~bolditalic(', Epithet, ')~bold(', Collection,')~bold(', Number,')')),
                color = 'black', parse = T) +
    geom_tiplab(aes(subset = (New == 0 & Type == 0),
                    label = paste0('italic(', Genus, ')~italic(', Epithet, ')', '~', Collection,'~', Number)),
                color = 'black', parse = T) +

    geom_text(aes(label= UFboot), vjust = 0, hjust = 1, size = 3) +
    xlim(NA,0.5)

```

Some taxa in the figure

# Answer
> 2 votes
I like both the methods and the study. `ggplot2` is R's standout package/library and the library `ggtree` looks great.

Anyway, I would know how to solve this in Python or Perl and I'll use the same logic for R. My suspicion is R is automatically type casting your metadata and therefore has no option but to remove the zero.

Simply convert your integer range into a string, viz.

```
 Number2 = as.character(Number)

```

or simply *in situ*

```
as.character(Number)

```

This doesn't make any difference to your output, who cares whether its and integer or string and forces R into retaining the zero, '0'.

---

Note below that tried and didn't work. A string should be immutable. The suspicion is the the '0' is being removed at an earlier point in the code, i.e. as soon as it type casts it - zero disappears.

What I would do some print outs of the variable between its import and writing to disk. What I suspect is the conversion to string needs to occur earlier.

---

**Second idea** The other possibility is to wrap the metadata in double/single quotes in the Excel file (or flat file). This way R cannot force it into an integer when it is imported and must retain it as a string. I'd be amazed if R could wriggle out of that.

---
Tags: phylogenetics, phylogeny, ggplot2, ggtree
---


# Question
Title: Why is the p-value significance threshold for HLA association tests $5*10^{-8}$?
Typically, a p-value of $p\<=5\*10^{-8}$ is used for **genome-wide** association testing, as there are roughly $1/p$ independent LD blocks in the human genome, so it correct for multiple testing across the genome.

HLA assocation tests, such as in this study, only look at a fraction of the genome, something like chr6:28–34 Mb, which is approximately 0.2% of the genome used to conduct a true 'genome-wide' association study. However, they still use the same genome-wide p-value threshold of $p\<=5\*10^{-8}$.

Is this not too conservative, since there is a much reduced multiple testing burden? Should the threshold not be something like $p/0.002$?

# Answer
Its a very good question and it is also a good study. What you are saying must be correct.

I thought it was possible they used the correct critical probability in the original submission and a reviewer requested it was changed. However I checked their medRxiv upload here and the authors clearly use this threshold prior their formal submission.

> We performed stepwise conditional analysis 479 to identify additional independent signals by adjusting for the most significant amino acid 480 position in each step until none met the significance threshold ( P = 5 10^-8)

Thus I don't know.

If this paper is central to your area of interest I would use the corrected critical probability threshold and re-assess and extend their conclusions because 5 decimal places is a huge a difference.
> 1 votes

# Answer
As a direct answer to your question, I would argue that the authors, in choosing this threshold, have demonstrated a lack of understanding of their process of analysis. But that answer is not a full answer, and it's a problematic answer because \[almost\] everyone who does population genetics does similar things. **This p-value is not appropriate, not merely because it is at a genome-wide significance level, but because *no* threshold is appropriate.**

To explain my perspective more, thresholding to identify "statistically-significant" results, in general, is not a good approach. Regardless of what threshold is used, p-values only indicate the confidence in a statistic's value (as it relates to fitting a particular statistical model), not importance or the size of the value.

The authors of this paper, as with many other authors who have written similar papers, make an obvious statistical error in the legend of figure 3, where they state: "One-field classical allele HLA-B\*57 (P = 9.84 × 10−138) (a) and amino acid position 97 in HLA-B (Pomnibus = 1.86 × 10−184) (b) showed the strongest association signal." This is an error, because they have only reported here the p-values, and p-values have no relevance to the strength or magnitude of a statistic.

Given this, it doesn't make sense to rank by p-value, and it especially doesn't make sense to display p-values on a graph and highlight p-value peaks as indicating association peaks. Those regions are *confidence* peaks.

In an attempt to provide some constructive light around this, I will answer the obvious follow-up questions are "How *should* p-values be used?" and "If p-values shouldn't be displayed or ranked, then what should?"

My preference is to use p-values as a loose filter, to exclude obviously incorrect results. In the case of a GWAS, I might, for example, choose to exclude any statistic that has a p-value above 0.1, leaving the remainder as *potentially* correct. It could possibly be argued that this is still thresholding. It is, in a sense, but importantly I am not saying that all values that pass the threshold are good, just that the ones that are excluded are all bad; I allow for some rubbish to end up in the good pile. I liken it to using a computer to find the things that look like diamonds, so that I don't waste my time inspecting dirt.

After that filtering is done, the p-values should be ignored, and the descriptive statistic used to rank, determine significance, association, or whatever the goal of the analysis is. For GWAS, a commonly-used value is beta, which can be positive or negative depending on the direction of association, and relates to the variant difference between groups.

If a list needs to be cut down, then this descriptive statistic should be used for thresholding, depending on the analysis. For example, it might be appropriate to display the 100 most extreme beta values, or threshold to exclude any beta values that have a magnitude less than 0.01, or display a manhattan-like plot with beta on the $y$ axis, and chromosome position on the $x$ axis.

There shouldn't be set rules, because every situation is different - there's a human element to this. What I'm saying is that it should be the job of the researcher, not the computer, to identify which data are significant. Make the data set small enough that it can be processed by a person, then let them do the rest. I think Hannah Fry explains this quite well in this talk:

> Now this is the stuff the algorithms are amazing at: looking for tiny, tiny clues in seemingly completely disconnected data sets.... Just incredibly, incredibly sensitive these algorithms, right? Just absolutely amazing what they can do. And yet remember that sensitivity does not make a perfect algorithm. That's only one half of the equation.... \[Humans\] are terrible at sensitivity, right? We are really, really bad. We miss things all the time. But, specificity, being specific, that's like our superpower.... And rather than choosing between human and machine, which is kind of the rhetoric that we get so much, why don't we exploit each other's strengths and just create much more of a partnership?.... Because the algorithm never gets tired, so let it trawl through all of that data and just highlight a few key areas of concern. And then the human never misdiagnoses. So they can come in, and just sweep up and have the final say.
> 1 votes

---
Tags: statistics, gwas, p-values
---


title: Why all values become 1 after dcast?

text:
# Question
Title: Why all values become 1 after dcast?
I am reshaping my dataframe as follows.

I want to set the `Wavelength column` as the column names with the code as follows.

```
dcast(unique(df1, by=cols), ID + Condition ~ Wavelength, value.var = 'mean')

```

But all the values `became 1s` as follows.

I would appreciate any help.

Yanpeng

# Answer
> 1 votes
That comment is right, it's a general programming question, and text would be much easier for people to help with rather than images.

But, since I ran into this not long ago...

Take a look at the message shown in blue just above the output table. My guess is you have some repeated ID/Condition/Wavelength combinations and it's defaulting to reporting the number of matching rows in the long-format data frame. You'll need to either make sure you have all unique combinations, or explicitly define what you expect it to do with multiple values going into a single output row/column position (like, median/mean/max/min/whatever). The unique() call won't do it since varying mean values would still make for unique input rows.

For example compare this:

```
> reshape2::dcast(data.frame(
+     Wavelength=c(664, 666, 662),
+     Condition="late",
+     ID=c(6221, 6222, 6223),
+     mean=c(106817, 110357, 112998)),
+   ID + Condition ~ Wavelength, value.var = "mean")
    ID Condition    662    664    666
1 6221      late     NA 106817     NA
2 6222      late     NA     NA 110357
3 6223      late 112998     NA     NA

```

with this:

```
> reshape2::dcast(data.frame(
+     Wavelength=c(664, 664, 666, 662),
+     Condition="late",
+     ID=c(6221, 6221, 6222, 6223),
+     mean=c(106800, 106817, 110357, 112998)),
+   ID + Condition ~ Wavelength, value.var = "mean")
Aggregation function missing: defaulting to length
    ID Condition 662 664 666
1 6221      late   0   2   0
2 6222      late   0   0   1
3 6223      late   1   0   0

```

---
Tags: r, statistics
---


# Question
Title: Running a list of apis via Python's request module
I am running the following Python code with `python requests` in my Jupiter notebook to get cancer oncology id's from EBI-EMBL.

```
import requests
cancer_all = requests.get('https://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0000376/children')
obi = cancer_all.json()

```

The above code is for one link. However I want to get results for 30 such links together.

Other examples of links:

```
https://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0000621/children
https://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0000637/children
https://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0000653/children
https://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0002402/children
https://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0002516/children
https://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0002813/children

```

**Questions**

1. Is there any way I can run all these api in one go?
2. Could I fetch the results (in json format) of these links in a single download. Thus there would be one json file for all the links?

# Answer
```
import requests
import json #??

apiString = "https://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0000621/children\nhttps://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0000637/children\nhttps://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0000653/children\nhttps://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0002402/children\nhttps://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0002516/children\nhttps://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0002813/children"

def requestIt():
    cancer_all = requests.get('https://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0000376/children')
    return cancer_all.json()

if __name__ == '__main__':
callList = split('\n')
obi = [requestIt(x) for x in callList]

# e.g.
print (obi[0]

# an alternative approach to list comprehension
# computationally faster but will not make much difference here
obi = list(map(requestIt, callList)

```

**Note** If you want to parallelise it, which might be what you meant by 'all at once', you instantiate `Pool()` from `multiprocessor` via `map` and use the same comprehension but this time the function `requestIt` will spread this over your available CPUs. Thats a little bit more complicated - namely getting the result back.

If the apis are in a file they will load into a list nicely:

```
with open ('apis.txt', 'r') as fin:
    callList = fin.readlines()

```

**However** you want to combine the json into a single json ... you `import json` then use `json.loads` to convert the json into a dictionary, combine the dictionaries and then use `json.dumps` to convert it back to json. The code for that is a separate question altogether and the comprehension might be better expressed as a dictionary comprehension (not sure).
> 1 votes

---
Tags: api, python, jupyter
---


# Question
Title: What does it mean by "topology" in the case of secondary and tertiary structures of proteins?
What does it mean by "topology" in the case of secondary and tertiary structures of proteins?

N.B. I am not talking about DNA/RNA/genomics. I am talking about protein folding, as it relates to protein modeling and simulation.

# Answer
> 2 votes
The word topology gets used in two different contexts.

## Common case: the geometric configuration of the protein (=fold)

In common non-computational biochemistry parlance it is just the standard English meaning of topology and just geometry, in the case of protein it means protein fold:

https://en.wikipedia.org/wiki/Protein\_fold\_class

Wiktionary (https://en.wiktionary.org/wiki/topology) say:

> The branch of mathematics dealing with those properties of a geometrical object (of arbitrary dimensionality) that are unchanged by continuous deformations (such as stretching, bending, etc., without tearing or gluing)

In the of protein, unchanged by small changes in length or amino acid sequence identity —no proteases or meat-glue ligase required!

## Computational (Bio)chemistry

Topology has a technical jargon definition in computational biochemistry and computational chemistry.

In computational chemistry how the connections of a molecule are configured, expressed as a graph network https://en.wikipedia.org/wiki/Graph\_theory, is referred to as topology. For example isopropanol and n-propanol have same number of heavy atom (*nodes*) and these will have the same hybridisation, but the network will have different connectivities (*vertices*). The topology definition generally runs off *atom types* more so than elements, these combine the element and hybridisation and other properties specific to a certain element in a molecule —e.g. amine vs. amide nitrogen vs. tryptophan ring nitrogens are different. The 3D embedding of said molecule is not part of its topology. Things get a bit problematic as ideal bond lengths, angles and dihedrals in dihedral space (i.e. not in cartesian space) go into the 'configuration' file (topology file) for several forcefield systems. And in others still there may be constraints or restraints specified.

In computational biochemistry, this holds true for ligands. But with proteins there is less. The connectivity within each amino acids is preset in stone: where one to use the allo isomer of threonine as opposed to the threo one would need to define a non-canonical amino acid. So by topology of a protein it is intended the connectivity of the amino acids in chains, e.g. primary sequence devoid of cartesian coordinates. The word topology appears often in the context of ligands: these require custom definitions (=topologies) required to use a ligand not parameterised in the force field used, a common source of trouble.

## See also

---
Tags: proteins
---