# Question
Title: Where does Snakemake send output that's generated by shell commands?
I'm wondering how Snakemake deals with output to `stdout` and `stderr` that is generated when running a rule with a `shell:` command.

I ran the following rule and only got informed that there was a non-zero exit code. However, I did not get a stack trace. To get access to the stack trace I had to run the shell command manually:

```
# Shortened for clarity
rule parse:
    input:
        sequences = rules.join_downloads.output
    output:
        sequences = "data/sequences.fasta"
    shell:
        """
        augur parse
        """

```

Why is that? And can I somehow make Snakemake output more detailed information when something goes wrong to help with debugging?

# Answer
> 1 votes
You usually define a log file like below, and redirect the stdout and/or stderr to it. Stdout is `>`, stderr is `2>` and both of them combined is `&>`.

```
rule parse:
    input:
        sequences = rules.join_downloads.output
    output:
        sequences = "data/sequences.fasta"
    log:
        my_output/my_log.log
    shell:
        """
        augur parse 2> {log}
        """

```

---
Tags: snakemake
---


title: Trinity assembly from many samples

text:
# Question
Title: Trinity assembly from many samples
When you combine samples for de-novo transcriptome assembly with Trinity, do you suggest limiting the number of reads for each sample? I had read one of Matthew MacManes' papers awhile back suggesting 40 million reads was a good number for Trinity assemblies of most metazoans.

However, I have 28 samples across 3 treatments that I want to use for DESeq and so, if I follow the advice at BioStars here, that's 28x40 million reads = 1.1 billion reads going into a Trinity assembler... This seems like a troublesome amount. How is it suggested to combine these samples in situations like this?

# Answer
Assuming these reads are from the same type of thing - and they should be if you're planning to use DESeq\[2\] on them - then it's fine to combine all the reads into a single Trinity assembly run.

All recent versions of Trinity will perform in silico normalization by default (based on kmer counts), so you shouldn't need to do a manual subsampling of reads. Read normalisation can be done separately; details on that can be found here:

https://github.com/trinityrnaseq/trinityrnaseq/wiki/Trinity-Insilico-Normalization

My vague memory is that the kmer counts are only used to get an idea of the total transcriptome size (i.e. for determining the average transcriptome coverage). If you want to get a more complete assembly, then consider only choosing samples that are most likely to have different overall transcript expression (e.g. different tissue types, treatments, etc.). If you've got 28 samples, you shouldn't need to worry about replicates; more variation is better.
> 1 votes

---
Tags: assembly, transcriptome
---


title: Converting a matrix array to matrix

text:
# Question
Title: Converting a matrix array to matrix
I have a matrix in R. When I run `class(my.matrix)` it returns:

```
"matrix" "array"

```

Why do I get "array" when in fact I want a plain matrix? The same `class` output is returned when running `as.matrix(my.matrix))`.

Why is that?

# Answer
> 5 votes
`class()` on a matrix will always return:

```
"matrix" "array"

```

Simple example:

```
> class(as.matrix(1))
[1] "matrix" "array"

```

This is because an array in R is a vector with multiple dimensions, and the matrix is a special type of array with exactly two dimensions, so every matrix is an array and this is what the `class` command on a matrix is telling us here, see http://adv-r.had.co.nz/Data-structures.html

# Answer
> 2 votes
From the return value of the `class` function, it looks like this already is a matrix. It just happens that it is considered by R to be both a matrix *and* an array.

---
Tags: r
---


# Question
Title: Output of allelic association doesn't write the rsID
I'm using plink (1.9b5) to do allelic association. My problem is that my output does not write the `rsID` (SNP), it just writes a dot.

Output:

```
CHR SNP  BP   A1  F_A      F_U      A2  P       OR
1   .    1014 A   0.18     0.07     C   0.22    3.03

```

My script:

```
plink --bed file.bed --bim file.bim --fam file.fam --make-bed --out output

```

Can someone please help me see what I'm doing wrong?

# Answer
> 0 votes
*Answer from @user438383, converted from comment:*

It seems likely your input .bim file probably doesn't have RSIDs either - is that the case?

---
Tags: snp, phylogenetics, ngs, plink
---


# Question
Title: GSEA: Computing an Enrichment Score (as described in the original paper)
The following paragraph comes from the GSEA paper that explains how to obtain an enrichment score for a gene set. Note that the enrichment score is computed **up to a position i** in L. There are 2 things I don’t understand about this:

1- first, the following image illustrating the method seems to suggest that that we score a gene set across the whole list of ranked genes, which conflicts with the idea of “up to a position i” in L

2- second, if we’re only considering genes up to a position i, then how are we going to capture the genes that are underexpressed (appear at the bottom of the ranked list) for phenotype A?

# Answer
Actually I believe I get it now. i is not a constant. We vary it from i=1 to N and consider the i for which $P\_{hit} - P\_{miss}$ is largest in magnitude
> 2 votes

---
Tags: gsea
---


# Question
Title: The confusion of using TPM (transcripts per million)
It is shown that TPM values are not suitable for DEG analysis but good for within-sample comparison since TPM normalized the gene length. My question is first: if TPM is not suitable for across sample comparison, then why are some heatmaps using TPM values? In what scenario are TPM values really useful? Should downstream analysis use TPM or raw count matrix?

When I analyzed sc-RNA seq data, the Seurat package uses normalized counts (which is CPM) for DEG, but TPM is not advised for differential testing. I'm not sure what I'm missing in the understanding. It seems very convoluted to me.

I just don't get the point that TPM is commonly used as an input for DEG testing by Seurat (Seurat findmarker function uses "data" slot, which is normalized counts), but TPM is not recommended for DEG testing for bulk.

# Answer
> 4 votes
Use raw counts as much as you can. Add in various relevant factors as covariates in DESeq2. RNAseq metrics have come a long way but are still misused by people because it's more convenient to compare some sort of normalized metric that give out relevant confounding factors, I guess.

# Answer
> 3 votes
The naive per-million scaling methods do not properly correct for the compositional bias between samples. This is especially true if the groups you compare are expected to be very different, e.g. different organs, see here demonstrated on some GTEx data: https://www.biostars.org/p/9465851/#9465854

Why do people use flawed metrics? Well, because some are very common. Per-million is easy to compute and does not change when new samples are added. The latter is convenient, and sometimes per-million might be good enough for visualization. I never do it though, I always use normalized (or vst) counts from DESeq2 or edgeR.

For differential analysis of bulk data one commonly uses raw counts which are then normalized internally by the established frameworks such as DESeq2, edgeR or limma-voom.

Single-cell data, if you consider each cell a replicate, often can go with simpler stats as the large sample such allows tests such as the T- or Wilcox test. Again, plain per-million scaling is probably bad, and alternative methods do exist, e.g. the deconvulution/sum-factor strategy in the `scran` package at Bioconductor. I cannot speak for Seurat, as I do not use it.

For most downstream applications such as visualization and clustering you probably want to use the properly normalized counts on the log2 scale. For something like heatmaps one would scale the data (Z-score) before plotting to show relative differences. I see no application for raw counts beyond feeding them into a DE testing framework that then normalizes them. The only applications that use raw countsbeyond the aforementioned DE testing frameworks that come to my mind are Combat-Seq (batch correction for bulk data) and sctransform (varaince stab/normalization/regression) for single-cell data.

# Answer
> 1 votes
For scRNA-seq, you don't want to normalize for gene length because the most popular 10X technology only sequences a the 3' or 5' end of the transcript. Hence CPM, which normalizes for sequencing depth alone, is popular for single cell data.

This is different from your average bulk RNA-seq, where you typically get full length transcripts (unless you're doing something like DRUG-seq).

---
Tags: rna-seq, scrnaseq, tpm
---


title: BEAUti not accepting taxa tip dates

text:
# Question
Title: BEAUti not accepting taxa tip dates
I'm using BEAUti as part of my pipeline and the sequence names include `|` to separate the metadata fields - with culture date being at the very end.

When I attempt to auto-configure the tip dates on the sequences however, some of the sequences will not update their date(raw value) and if I double-click and manually change the date it will then reverse back to 0.

Any idea what may be preventing me from manually correcting the dates or auto-config to ignore the culture date after the last `|`?

# Answer
> 1 votes
I found the solution, one of my Seq IDs had a `,` in one of the meta-data fields, this appears to have thrown off the software and all subsequent Seq IDs could not be auto-configured.

Simply removing the `,` solved the issue.

---
Tags: multiple-sequence-alignment, beast
---


title: Remove from Multi-FASTA by Sequence ID

text:
# Question
Title: Remove from Multi-FASTA by Sequence ID
I want to remove sequence VRE32514 – it doesn’t belong and thus is the reason it lacks additional metadata. However I tried implementing this code from a similar question:

```
awk 'BEGIN{while((getline<">VRE32514")>0)l[">"$1]=1}/^>/{f=!l[$1]}f' *.fasta | grep -v ">VRE32514" *.fasta

```

With no success and the proposed solution in R using `seqinr` doesn't seem to work either for myself.

```
#load package
library(seqinr)
#load file containing sequences
data<-read.fasta("test.fasta")
#create a vector containing species names: these are the last part of the string
vec.names<-unlist(lapply(strsplit(names(data), "|"), function(x)x[length(x)]))
#find names to keep: indices which are not in the species to remove
species.to.remove<-c(">VRE32514")
vec.tokeep<-which(! vec.names %in%  species.to.remove)
#write the final output
write.fasta(sequences=data[vec.tokeep], names=names(data)[vec.tokeep], file.out="output.fasta")

```

Below is an example of what my Multi-FASTA appears as and the desired output:

```
>VRE32491|736|PUH-10C|Blood|2016-12-07
ATGAGATCAGAAAAAGAAATGATGGATTTAGTACTTTCTTTAGCAGAACAGGATGAACGT
ATTCGAATTGTGACCCTTGAGGGGTCACGCGCAAATATTAATATACCTAAAGATGAATTT
>VRE32493|1471|PUH-10N|Tissue/Surgical|2016-12-08
CAGGATTATGATATTACATATTTTGTAAGTGATATAGAACCGTTTATATCTAATGATGAC
TGGCTTAATCAATTTGGGAATATAATAATGATGCAAAAGCCGGAGGATATGGAATTATTC
CCACCTGAAGAAAAGGGATTTTCCTATCTTATGCTATTTGATGATTACAATAAAATTGAT
>VRE32503|1471|PUH-11N|Wound|2017-01-05
CTTACCTTATTGCCCTTGGAAGAGTTAGATAATTACCTAAAGGGCGATAAATTAATAAAG
GTTCTAATTGATAAAGATTGTAGAATTAAAAGGGACATAGTTCCGACTGATATAGATTAT
>VRE32504|1471|PUH-EMEP|Blood|2017-01-10
CATGTAAGAAAGCCAAGCGCAAGGGAGTATGATGATTGCTGCAATGAATTTTGGAATGTA
ACACCTTATGTTATTAAAGGATTGTGCCGTAAGGAAATTTTATTTGCTATTGATCATTTT
>VRE32514
AATCAGATTGTTCGCCATGAGCTGCTGAGAATGATATCATGGAAGGGCGGCATCGAAACA
GGCTTTAAATTAAGTGTAGGCAAGAACTATAAGTTTATTGAAAGGTATATATCCGAGGAT

```

Desired Output:

```
>VRE32491|736|PUH-10C|Blood|2016-12-07
ATGAGATCAGAAAAAGAAATGATGGATTTAGTACTTTCTTTAGCAGAACAGGATGAACGT
ATTCGAATTGTGACCCTTGAGGGGTCACGCGCAAATATTAATATACCTAAAGATGAATTT
>VRE32493|1471|PUH-10N|Tissue/Surgical|2016-12-08
CAGGATTATGATATTACATATTTTGTAAGTGATATAGAACCGTTTATATCTAATGATGAC
TGGCTTAATCAATTTGGGAATATAATAATGATGCAAAAGCCGGAGGATATGGAATTATTC
CCACCTGAAGAAAAGGGATTTTCCTATCTTATGCTATTTGATGATTACAATAAAATTGAT
>VRE32503|1471|PUH-11N|Wound|2017-01-05
CTTACCTTATTGCCCTTGGAAGAGTTAGATAATTACCTAAAGGGCGATAAATTAATAAAG
GTTCTAATTGATAAAGATTGTAGAATTAAAAGGGACATAGTTCCGACTGATATAGATTAT
>VRE32504|1471|PUH-EMEP|Blood|2017-01-10
CATGTAAGAAAGCCAAGCGCAAGGGAGTATGATGATTGCTGCAATGAATTTTGGAATGTA
ACACCTTATGTTATTAAAGGATTGTGCCGTAAGGAAATTTTATTTGCTATTGATCATTTT

```

# Answer
> 2 votes
You can use `samtools faidx` to first index your file and then retrieve only those sequences you specify. For this `grep` all the sequence headers you want and then remove the unwanted ones, followed by retrieval of the fasta entries with `faidx`:

```
#/ index:
samtools faidx your.fa

#/ run:
grep '^>' test.fa \             # grep names
| grep -v '>VRE32514' \         # remove unwanted
| tr -d '>' \                   # remove ">" from name
| xargs samtools faidx your.fa  # get the entries you want

```

Output:

```
>VRE32491|736|PUH-10C|Blood|2016-12-07
ATGAGATCAGAAAAAGAAATGATGGATTTAGTACTTTCTTTAGCAGAACAGGATGAACGT
ATTCGAATTGTGACCCTTGAGGGGTCACGCGCAAATATTAATATACCTAAAGATGAATTT
>VRE32493|1471|PUH-10N|Tissue/Surgical|2016-12-08
CAGGATTATGATATTACATATTTTGTAAGTGATATAGAACCGTTTATATCTAATGATGAC
TGGCTTAATCAATTTGGGAATATAATAATGATGCAAAAGCCGGAGGATATGGAATTATTC
CCACCTGAAGAAAAGGGATTTTCCTATCTTATGCTATTTGATGATTACAATAAAATTGAT
>VRE32503|1471|PUH-11N|Wound|2017-01-05
CTTACCTTATTGCCCTTGGAAGAGTTAGATAATTACCTAAAGGGCGATAAATTAATAAAG
GTTCTAATTGATAAAGATTGTAGAATTAAAAGGGACATAGTTCCGACTGATATAGATTAT
>VRE32504|1471|PUH-EMEP|Blood|2017-01-10
CATGTAAGAAAGCCAAGCGCAAGGGAGTATGATGATTGCTGCAATGAATTTTGGAATGTA
ACACCTTATGTTATTAAAGGATTGTGCCGTAAGGAAATTTTATTTGCTATTGATCATTTT

```

In `R` one could use `Biostrings`:

```
library(Biostrings)

fasta <- readDNAStringSet("path/to/your.fa")
fasta_new <- fasta[!names(fasta) %in% c("VRE32514")]

#/ save to disk:
writeXStringSet(fasta_new, "path/to/outputfile.fa")

```

# Answer
> 0 votes
Here's a solution using Python:

from Bio import SeqIO

sequences = SeqIO.parse('test.fasta', 'fasta')
filtered = [seq for seq in sequences if 'VRE32514' not in seq.id]

with open('output.fasta', 'wt') as output:
    SeqIO.write(filtered, output, 'fasta')
```

---
Tags: fasta, linux, shell, awk, multiple-sequence-alignment
---


# Question
Title: A unified database for CNV, SNP, Indel and MSI
I am looking for a database or different databases where I can find information on different gene variants for a gene. As an example if I enter PPARG, I could be able to see SNPs, CNV, InDels and MSI.

In my understanding the information on SNPs can be downloaded from dbSNPs, MSI from ClinVar, CNV & InDels from dbVar. Please correct me if I am wrong.

Is there a database that houses all of this information at under one roof?

Your suggestions are very welcome.

Thank you

# Answer
To Answer my own question, one can use variation viewer to get all the information :)

DisGeNet is another resource.
> 1 votes

---
Tags: snp, mutations, cnv, indel
---


# Question
Title: Cell Ranger) Error during performing a test run
I am interested in the analysis of single-cell RNA-sequencing data. I installed Cell Ranger according to the instruction (https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/tutorial\_in#download). I know, it should be easy, but I found an error occurring after test running as below.

```
[error] Pipestance failed. Error log at:
tiny/SC_RNA_COUNTER_CS/SC_RNA_COUNTER/_BASIC_SC_RNA_COUNTER/COUNT_GENES/fork0/chnk0-u6c389ca928/_errors

```

Within that file, it shows the following error:

```
Job failed in stage code. signal: killed.

```

So, what do you think about this problem?

Additional information:

* a checksum test (to check the integrity of the downloaded files) was successful
* I performed a test run on a personal laptop, mac.
* I used SSH command to connect the server, but we are sharing the server with our lab members.

I found an error as follows:

```
WARNING: configured to use 112GB of local memory, but only 80.5GB is currently available.

```

However, at this time with unknown reasons, the installation is successfully completed, in spite of the warning. Do you have idea?

Thank you in advance.

# Answer
*Answer from PPK, converted from comments:*

The memory warning is only letting you know that the upper limit you set on memory (112 gb) was higher than the available memory. Do you share the server with other people? because the minimum amount of memory for Cell Ranger is 64 GB. In your initial run there may have been less than 64 GB free.
> 1 votes

---
Tags: scrnaseq, software-installation, cellranger
---


# Question
Title: General question: How to assign my SNP (Wheat, Illumina 15K) to specific chromosomes?
I have 2000 lines of wheat genotyped with a 15K illumina chip. I also have the name of the SNPs. I would like to know how could I find out which SNPs are in each chromosomes.

# Answer
> 0 votes
*Answer from @devon-ryan, converted from comments:*

Look in Biomart to see if there's SNP rsids for wheat.

---
Tags: snp, chromosomes
---


title: Phasing partially phased genomic data

text:
# Question
Title: Phasing partially phased genomic data
I am trying to create an integrated callset with coding regions from phased exomes and noncoding regions from a genotyping array. I would like to merge the data together in the correct phase. Since the exomes are already phased with read informed phasing, I don't want to rephase them because some of the phasing changes after doing this. Is there a way to statistically phase genetic data that has already been partially phased?

My data is in VCF format.

Thanks

# Answer
*Answer from @user438383, converted from comment:*

I would recommend checking out shapeit4 as it is able to integrate statistical phasing and read-based phasing like you are after. The manual page for WhatsHap might also be useful for you:

https://whatshap.readthedocs.io/en/latest/index.html
> 1 votes

---
Tags: genomics, genotyping, exome, phasing
---


title: Pol_II genome wide occupancy analysis

text:
# Question
Title: Pol_II genome wide occupancy analysis
I would like to know how to perform this global Pol II occupancy which is shown here:

Paper

Figure 1G

So in this figure they are analysing the Pol II occupancy in Control,PAF1 KD and LEO1.

I have similar data.

> In order to determine whether PAF1 has a direct effect on the regulation of Pol II at the promoter-proximal region, ChIP-seq was performed. Heatmaps ranked by decreasing Pol II occupancy show a positive correlation with PAF1 and Pol II (Figure 1G). To corroborate these findings, we also performed ChIP-seq with LEO1, as this was the only other antibody that we tested that gave any significant ChIP-seq signal. Like PAF1, LEO1 is found throughout gene bodies with particular enrichment around the TSS (Figure 1G). Interestingly, LEO1 was not broadly lost from chromatin upon PAF1 knockdown, but the distribution was altered, including reduced occupancy in the promoter-proximal region

So what I have done so far in terms of analysis is I did peak-calling using homer for my control and KD condition along with their Pol II chip.I was going for this differential binding analysis but this is now how i Pol II data is done after reading the paper I got. In homer there is this function analyzeRepeats.pl which is provided in this section **Measuring Gene Expression in Exons vs. Gene Bodies.** telling genes *(default) - Counts tags on the full gene body (TSS to TTS).This is useful for GRO-Seq where we expect coverage across the entire transcript.Can also be used to quantify H3K36me3 or PolII ChIP-Seq.*

How do I proceed and what tool should I use to show the Pol II occupancy at different condition and as of now I don't have respective rna seq data, do i need expression data of genes to perform occupancy analysis.

# Answer
*comment discussion between @krushnach-chandra and @stupidwolf (in chronical order):*

**SW:** No you don't need gene expression data.. You can just count the reads from polII chip (normalized by some method of your choice) and perform the analysis

**KC:** are you suggesting I can compare the Pol II chip between control and Knock out with their respect IP data?

**SW:** more or less yes...if you are interested in genomewide effects.. For example, if the treatment causes the polymerase to stall, you will see more reads in the promoter compared to the gene body

**KC:** so naive question do i need those Input and IgG at all for the comparison?

**SW:** good question.. depends on how you analyze it. If you always normalize within the gene, for example, you calculate a coverage normalized to that in the gene, or some meta gene thing, you don't need.

**SW:** If you do a differential gene expression based analysis, in theory you don't really need that, because you will always compare between your two conditions right. The IgG is good for checking you don't have any odd stuff... You can see whether there are genes where the control is enriched over your samples

**KC:** "calculate a coverage normalized to that in the gene, or some meta gene thing" any example can you cite me..so that i can follow and proceed ,as it looks like binding data are really complicated compared to rna seq

**KC:** "If you do a differential gene expression based analysis," you saying i can take out counts and do diff analysis for chip seq data using deseq2? but how do i get the information regarding if Pol II activity is going or stalled ? and how do i find the occupancy from the differential result

**SW:** welcome to the world of epigenomics!!! maybe something like this

**SW:** As for differential gene expression, yes. If you do this separately for promoters and gene body, you can see that the fold change will be different for a gene's promoter vs gene body, in situations where there's stalling, or reduced elongation etc.. I think it's a good way to account for the variation

**KC:** your insights are good start for me to run the analysis

**KC:** mundane question how do i get promoter and gene body information is it after annotation which should be done after running analysis ?

**SW:** Meaning? I normally align the reads, then for TSS and gene body you have to define it yourself using the transcript or gene annotation. Normally like +/- 200 bp around TSS for promoter and gene body is rest of it?

**KC:** "Normally like +/- 200 bp around TSS for promoter and gene body is rest of it" so i can do that in both chipseeker as well as homer?
> 1 votes

---
Tags: chip-seq, deeptools, homer
---


title: Finding difference between groups

text:
# Question
Title: Finding difference between groups
I have mutational catalogues of 4 samples like below

```
> head(out)
        Responders1 Nonresponders1 Responders2 Nonresponders2
A[C>A]A       11744          13546        2897           3655
A[C>A]C        5144           7172        2295           2163
A[C>A]G         939           1257         290            279
A[C>A]T        6065           7997        2078           2088
C[C>A]A        8969          11155        2055           2582
C[C>A]C        4050           6657        1367           1173
>

 > dim(mut)
[1] 96  4
>

```

*Note: A\[C\>A\]A - Single nucleotide base substitution*

I want to look if mutational substitutions are different between samples

I have tried

```
wilcox.test(mut$Responders2,mut$Nonresponders2)

t.test(mut$Responders2,mut$Nonresponders2)

var.test(mut$Responders2,mut$Nonresponders2)
fisher.test(mut$Responders2,mut$Nonresponders2,simulate.p.value=TRUE)

kruskal.test(mut$Responders2,mut$Nonresponders2)

```

To find a difference but I think you may know a better test to find any difference here

This is base shift (mutation, substitution across genome of each patient). I guess mutations happen randomly but I don't know if mutational caller algorithm also pick them by randomness. I expect the shift of A (Adenine) to C (Cytosine) in a group of these patients be higher than the other group. This higher rate likely should be due to some background than randomness picked by algorithm. So I think null hypothesis is that rate of base shifting is not the same in both group. I also don't expect a normalised distribution of base shifting rather a poison distribution.

Can you help please?

Thanks

# Answer
*Answer from @fabio-marroni, converted from comment:*

You might try to use mutationalPatterns and see if the responders cluster together and separated from nonresponders.
> 1 votes

---
Tags: r, wgs, mutations
---


# Question
Title: Heatmap annotation with value and colour gradient
I am trying to use ComplexHeatmap to plot heatmap with annotation. The code that I have tried

```
library(ComplexHeatmap)
library(circlize)

set.seed(123)
mat = matrix(rnorm(25), 5)
rownames(mat) = paste0("R", 1:5)
colnames(mat) = paste0("C", 1:5)

col_fun = colorRamp2(c(0, 2, 5), c("white", "blue", "red"))
data_foo = c(0.1, 3, 1.6, 2.8, 4)

ha = HeatmapAnnotation(which = "row", foo = data_foo, col = list(foo = col_fun))
Heatmap(mat, name = "mat", right_annotation = ha)

ha1 = HeatmapAnnotation(which = "row", foo = anno_text(data_foo))
Heatmap(mat, name = "mat", right_annotation = ha1)

```

Can I obtain the annotation with colors varying as per the annotation values and also display the values in the annotation? i.e is it possible to combine the annotation in both figures?

# Answer
*Answer from @haci, converted from comment:*

You can use the `pch` parameter in `anno_simple()` but you can only use a single character or a digit. It is not limited to symbols only, you can provide a vector to `pch`, but only the first character would be printed. If by chance your values are \[0-9\] you can use `pch` as kind of hack for your purposes above.
> 0 votes

---
Tags: r, heatmap
---


# Question
Title: Arm-level SCNA does not change mRNA expression level?
It is almost common sense that local amplification of gene increase the its mRNA expression . I recently heard that large scale SCNA (such as arm-level amplification) does not increase the mRNA expression level of gene within the region. I did some analysis (SCNA is determined by SNP array and NGS, mRNA expression is measured by RNA-Seq) and found that this is indeed the case. I further did some research and did not found any explanations. it seems like a common acknowledgement. Is there any technical bias or biological explanation that I am missing here? Or is it something that has not have a answer yet?

I have no clue about this at all. Can you recommend this papers on this?

# Answer
*Answer from @devon-ryan, converted from comments:*

The common presumption in such cases is that there's a feedback loop on most genes to control their expression level. It doesn't have to be specific to arm-level expression.
> 1 votes

---
Tags: rna, cnv
---


# Question
Title: Why the gene symbol aliases have significantly different sequences?
"Aliases" or "Synonyms" should represent the same gene but with different names. But when I try to find the sequence of one alias and try to match with another alias they are significantly different.

During bioinformatics analysis, I need to convert NAT5 alias ID to its gene symbol ID (NAA20, NAA50). When I tried to match the sequences of NAA20 and NAA50, they are very different but still NAT5 is an alias of both NAA20 and NAA50.

I have to convert one id to another id (ENSEMBL). These gene symbol and alias is so confusing.

I am looking for NAT5 and when I search for NAT5 in ENSEMBL, it shows both NAA20 and NAA50.

From UniProt the canonical sequence of both the proteins:

https://www.uniprot.org/uniprot/P61599#sequences

https://www.uniprot.org/uniprot/Q9GZZ1#sequences

I am confused which one to chose?

# Answer
> 5 votes
At Ensembl, we categorise synonyms as anything that a gene might also be known as. This includes older names for them, since those names will be in the literature, including where a gene has been split in two.

# Answer
> 2 votes
*Answer from @devon-ryan, converted from comment:*

Both, that's why NAT5 is no longer a human gene symbol.

NAT5 does not equal NAA20, it's an out-dated name for it. Names change over time as people realize that there are more genes for something than originally thought.

---
Tags: ensembl, identifiers, entrez
---


title: working with cutadapt

text:
# Question
Title: working with cutadapt
I'm working with ion torrent data where I apply the program cutadapt. I'm analyzing ITS seq data, as well as Matk. When I'm using cutadapt, I search the information for this genes, and built new files accordingly to his information. For instance, when I search for ITS, I generate a file about this gene, and the same applies to the Matk gene. However, more recently I randomly picked a fasta sequence from ITS, blast it and received the information about the matk gene. Is this supposed to occur? Thank you for the help and so sorry for the naive question.

# Answer
*Answer from @m, converted from comment:*

It's definitely not supposed to happen. There appears to be a deep rooted classification bug, where matk is misclassified as ITS. Unpicking this will be hard.
> 0 votes

---
Tags: python, sequence-analysis, cutadapt
---


# Question
Title: Can the header line contain arbitrary text in FASTA format?
The following is a valid .fasta file content:

`>HSBGPG Human gene for bone gla protein (BGP) GGCAGATTCCCCCTAGACCCGCCCGCACCATGGTCAGGCATGCCCCTCCTCATCGCTGGGCACAGCCCAGAGGGT ATAAACAGTGCTGGAGGCTGGCGGGGCAGGCCAGCTGAGTCCTGAGCAGCAGCCCAGCGCAGCCACCGAGACACC ATGAG.`

Is this also?

`>Arbirary_Name_iJustCameUP_with_and_other_local_identifiers GGCAGATTCCCCCTAGACCCGCCCGCACCATGGTCAGGCATGCCCCTCCTCATCGCTGGGCACAGCCCAGAGGGT ATAAACAGTGCTGGAGGCTGGCGGGGCAGGCCAGCTGAGTCCTGAGCAGCAGCCCAGCGCAGCCACCGAGACACC ATGAG.`

# Answer
NCBI's FASTA format description:

https://www.ncbi.nlm.nih.gov/genbank/fastaformat/

NCBI's BLAST page describing valid FASTA input:

https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE\_TYPE=BlastDocs&DOC\_TYPE=BlastHelp

NCBI's SNP page describing FASTA format:

https://www.ncbi.nlm.nih.gov/projects/SNP/snp\_legend.cgi?legend=fasta

UniProt's FASTA-header description:

https://www.uniprot.org/help/fasta-headers

Wikipedia reference (indicating format origin):

https://en.wikipedia.org/wiki/FASTA\_format

FASTA program (origin of format, per Wikipedia):

https://fasta.bioch.virginia.edu/wrp\_fasta/fasta\_guide.pdf

Harvard *PolyPhen* page describing FASTA format:

http://genetics.bwh.harvard.edu/pph/FASTA.html

FASTA file format (from Pacific Biosciences):

https://pacbiofileformats.readthedocs.io/en/3.0/FASTA.html
> 2 votes

# Answer
I am not aware of an official fasta format description. The only constraint that I know of is that there should not be any whitespace after `>`. Apart from that particular position, whitespaces are allowed in the header line.

So, regarding your particular example, both are valid fasta sequences. When in doubt, you can use `SeqIO` from `Biopython`, if you can parse your file with the following code, it  is should be a valid fasta file.

```
from Bio import SeqIO

with open("example.fasta") as handle:
    for record in SeqIO.parse(handle, "fasta"):
        print(record.id)

```

Edit per @Chris\_Rands' comment

The code below does the same as above, meaning that `SeqIO.parse()` takes care of opening and closing the file.

```
from Bio import SeqIO

for record in SeqIO.parse("example.fasta", "fasta"):
    print(record.id)

```
> 4 votes

---
Tags: fasta
---


# Question
Title: Seurat Dimplot with different clustering IDs
In Seurat metadata I have assigned cells to some cell types with different resolutions

I have added the cluster identities to the object via `Idents`:

```
 Idents(control) <- a$predicted.celltype.l1

```

I then added another clustering in the same way:

```
 Idents(control) <- a$predicted.celltype.l2

```

Dimplot will now show the new identities:

How I can ask Dimplot to plot the clustering I want for example if I want to show `predicted.celltype.l1` again (something like changing the default)?

# Answer
The thing to tweak is the `group.by` argument of the `Dimplot` function.

It looks like you've already tried this in your attempts at solving a different problem.

The main thing you need to be aware of is that each new clustering needs to be saved as a separate metadata variable in order to do this. As I mentioned in one of my previous answers, I like to use the `[[X]]` notation for this:

```
control[["predicted.celltype.l1"]] <- a$predicted.celltype.l1
control[["predicted.celltype.l2"]] <- a$predicted.celltype.l2

```

Then you should be able to use both in the `group.by` argument:

```
DimPlot(control, reduction = "umap", group.by = "predicted.celltype.l1")
DimPlot(control, reduction = "umap", group.by = "predicted.celltype.l2")

```
> 1 votes

---
Tags: scrnaseq, seurat, 10x-genomics
---


title: What does this Seurat argument mean

text:
# Question
Title: What does this Seurat argument mean
I have extensively read about `percent mito` in Seurat but I got more and more confused

Let's say we want to keep cells with `maximum 20%` of their expressed genes are mitochondrial genes, so shall I set `percent.mito < 20` or `percent.mito < 0.2` or `percent.mito < 80` ?

# Answer
> 2 votes
Checkout the first Seurat tutorial. To filter cells with \>20% mitochondrial counts:

```
obj[["percent.mito"]] <- PercentageFeatureSet(obj, pattern = "^MT-")
obj <- subset(obj, subset = percent.mito < 20)

```

---
Tags: scrnaseq, seurat, 10x-genomics
---


title: Filtering cells from values in metadata

text:
# Question
Title: Filtering cells from values in metadata
I metadata slot of Seurat object I have mapping score of each cell to a reference PBMC data like

```
pbmc <- subset(object, subset = object@meta.data$mapping.score > 0.5)
Error in FetchData(object = object, vars = unique(x = expr.char[vars.use]),  :
  None of the requested variables were found:
>
>
>

```

I want to only keep cells with mapping score \> 0.5

Any help please?

# Answer
> 1 votes
```
your_new_object <- subset(your_object, subset = mapping.score > 0.5)

```

---
Tags: scrnaseq, seurat, 10x-genomics
---


# Question
Title: Convert paired-end BAM into a single-end BAM and keep all the reads
I want to call peaks using MACS2 on my paired-end ATAC-seq BAM and I am interested in cutting sites. So I make use of the -BAM option (opposed to the -BAMPE) option. However this throws away all the mates, and all the reads of which the mate did not align (that's more than 50%!!!). My idea was to align paired-end, and then remove all the information from the BAM that it is paired-end. So my question is how do I:

**Convert paired-end BAM into a single-end BAM and keep all the reads**

I made a start based on the SAM-format specification, but when I run the result of this script through markduplicates I get the error:

```
ERROR: Record 1, Read name SRR6245286.60336908., MRNM should not be set for unpaired read.

```

I can not find what this MRNM is...

This is my code I used to generate the faulty *'single-end'* bam, but I am open to suggestions (e.g. bash, preferably no R).

```
import pysam

paired =         1
proper_pair =    2
mate_unmapped =  8
mate_reverse =   32
first_in_pair =  64
second_in_pair = 128

bam_in = pysam.AlignmentFile("/vol/atacpipeline/bwa/GSM2837489-danRer7.samtools-coordinate.bam", "rb")
bam_out = pysam.AlignmentFile("/vol/atacpipeline/bwa/GSM2837489-danRer7_SE.samtools-coordinate.bam", "wb", template=bam_in)
for line in bam_in:
    if line.flag & second_in_pair:
        line.pos += line.template_length
    line.next_reference_id = 0
    line.next_reference_start = 0
    line.flag &= ~(paired + proper_pair + mate_unmapped + mate_reverse + first_in_pair + second_in_pair)

    bam_out.write(line)

```

# Answer
> 2 votes
Like Devon pointed out, most likely you should sort out whether the files have been marked for duplicates correctly.

You can also use samtools rmdup too

```
samtools rmdup <bamfile> <output>

```

Also, I find it very odd that 50% of the reads are lost when you keep only properly paired reads.

I would suggest removing the duplicates and then converting the bam file to a bed file. Macs2 can take in the bed format. Any extremely crude python script that converts a bam into a bed:

```
import pysam
samfile = pysam.AlignmentFile("test.bam", "rb")
open("test.bed", 'w').close()
f = open("test.bed", "a")
STRAND = ["+","-"]
for read in samfile.fetch():
    STR = STRAND[int(read.is_reverse)]
    BED =[read.reference_name,read.pos,read.reference_end,".",read.mapq,STR]
    f.write('\t'.join([str(i) for i in BED]))
    f.write('\n')
f.close()

```

You can also try the conversion with bedtools or definitely refine that script above with pandas etc.

# Answer
> 2 votes
MRNM stands for "Mate reference index". So Picard found something in the RNEXT field which should be set only for paired-end reads but the rest of the file looks like single-end.

The problematic line in your code is:

```
line.next_reference_id = 0

```

This sets the RNEXT SAM field to whatever Pysam stores as a reference with index 0 (`next_reference_id`). This is **not** the same as reference name (`next_reference_name`). In single-end SAM file this field should be set to `*`. You cannot set `*` with `next_reference_id` but you can set it with `next_reference_name`:

```
line.next_reference_name = "*"

```

I would also recommend to change `next_reference_start` line to:

```
line.next_reference_start = -1

```

just to be sure Pysam handles this field as not set (I tried to put there `0` Pysam changed it to `1` for some reason). I would also add line:

```
line.template_length = -1

```

because this field should also be empty in single-end SAM.

This worked in my case. If Picard keeps complaining you can add `VALIDATION_STRINGENCY=LENIENT` or `SILENT` to your Picard command to suppress exiting of the command at errors.

---
Tags: pysam, macs2
---


title: Merging two dataframes in R

text:
# Question
Title: Merging two dataframes in R
I have a file with cluster number of a scRNA-seq and corresponding annotated cell type like

```
> dput(head(anno))
structure(list(cluster = 1:6, celltype = c("T:CD4+NAIVE", "T:CD4+NAIVE",
"T:CD8+NAIVE", "NK:CD56+16+3+NKT", "T:CD8+NAIVE", "B:")), row.names = c(NA,
6L), class = "data.frame")

> unique(anno$cluster)
 [1]  1  2  3  4  5  6  8  9 11 12 15 16 19 20 21 22 23 24 25 27 28 29 30 31 35 36 37 38 39 40 41 42 43 44 45 46
[37] 47 48 49 50 51 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 80 81 82 83 84
[73] 85 86 87 88 89 90 91 92 93 94 95
> unique(anno$celltype)
 [1] "T:CD4+NAIVE"        "T:CD8+NAIVE"        "NK:CD56+16+3+NKT"   "B:"                 "NKT"
 [6] "T:CD8+EM"           "T:CD4+CM"           "mDC:"               "MACROPHAGE:"        "MONOCYTE:precursor"
[11] "NEUTROPHIL:"        "NK:CD56+16+3-"      "MAST:"              "MONOCYTE:"          "pDC:"
[16] "T:Reg"

```

and another data frame with cell IDs and corresponding cluster number of each cell

```
> dput(head(df2))
structure(list(cluster = c(1L, 1L, 1L, 1L, 1L, 1L), cellid = c("cell17203",
"cell17205", "cell17206", "cell17207", "cell17208", "cell17209"
)), row.names = 15498:15503, class = "data.frame")

> unique(df2$cluster)
 [1]  1  2  3  4  5  6  8 12 16 19 21 23 27 29 30 35 36 40 53 55 92 11 24 38 41 46 90  9 15 20 47 74 57 39 42 25
[37] 56 64
>
> dim(df2)
[1] 16056     2
>

```

I want to create a cell type column in the second data frame (df2) in which relate cell type from the first data frame to cluster number for instance `cluster 1` is `T:CD4+NAIVE`

I have tried these

```
> df2[anno, on = "cluster"]
Error in `[.data.frame`(df2, anno, on = "cluster") :
  unused argument (on = "cluster")
> df2[c, on = anno$cluster]
Error in `[.data.frame`(df2, anno, on = c$cluster) :
  unused argument (on = anno$cluster)
> df2[c, on = df2$cluster]
Error in `[.data.frame`(df2, anno, on = df2$cluster) :
  unused argument (on = df2$cluster)

```

# Answer
You used `data.table` notation without actually having data.table:

```
library(data.table)
setDT(anno)
setDT(df2)
df2[anno, on = "cluster"]

```

```
    cluster    cellid         celltype
 1:       1 cell17203      T:CD4+NAIVE
 2:       1 cell17205      T:CD4+NAIVE
 3:       1 cell17206      T:CD4+NAIVE
 4:       1 cell17207      T:CD4+NAIVE
 5:       1 cell17208      T:CD4+NAIVE
 6:       1 cell17209      T:CD4+NAIVE
 7:       2      <NA>      T:CD4+NAIVE
 8:       3      <NA>      T:CD8+NAIVE
 9:       4      <NA> NK:CD56+16+3+NKT
10:       5      <NA>      T:CD8+NAIVE
11:       6      <NA>               B:

```
> 1 votes

# Answer
For merging dataframes, I find it easiest to use the tidyverse / dplyr functions inner/full/left/right\_join. See the "Data Transformation Cheatsheet" on this page. For the merge that @user438383 has mentioned, this would be a *left\_join*:

```
library(tidyverse)
anno <- data.frame(cluster = 1:6, celltype = c("T:CD4+NAIVE", "T:CD4+NAIVE",
                   "T:CD8+NAIVE", "NK:CD56+16+3+NKT", "T:CD8+NAIVE", "B:"))

df2 <- data.frame(cluster = c(1L, 1L, 1L, 1L, 1L, 1L),
                  cellid = c("cell17203", "cell17205", "cell17206",
                             "cell17207", "cell17208", "cell17209"),
                  row.names = 15498:15503)

> left_join(anno, df2)
Joining, by = "cluster"
   cluster         celltype    cellid
1        1      T:CD4+NAIVE cell17203
2        1      T:CD4+NAIVE cell17205
3        1      T:CD4+NAIVE cell17206
4        1      T:CD4+NAIVE cell17207
5        1      T:CD4+NAIVE cell17208
6        1      T:CD4+NAIVE cell17209
7        2      T:CD4+NAIVE      <NA>
8        3      T:CD8+NAIVE      <NA>
9        4 NK:CD56+16+3+NKT      <NA>
10       5      T:CD8+NAIVE      <NA>
11       6               B:      <NA>

```
> 5 votes

---
Tags: r, scrnaseq, seurat, 10x-genomics
---


# Question
Title: Low pass sequencing has been reported to detect common variants. How low can one go and get reliable data? Is 2X pass sequencing analysis possible?
I would like to use low pass sequencing to replace a genotyping chip to be able to detect variants up to 0.1 % allele frequency in available population data. What is the minimum depth I can opt for to retain the confidence of base calls? Can GATK be used to analyse this data?

# Answer
> 2 votes
This is a good idea - low pass sequencing has shown to be better than genotype arrays for things such as GWAS or QTL mapping.

The primary way of processing low-pass sequencing data is to use imputation - partciularly methods such as GLIMPSE or QUILT and a large reference panel such as the HRC or TopMed in order to refine the genotype calls.

This is a figure from the GLIMPSE paper showing how you can retain variants of different frequencies at different levels of sequencing depth.

And a similar image from the QUILT paper:

You might also want to check out some of the methodology in this paper http://m.genome.cshlp.org/content/31/4/529

---
Tags: ngs, gwas, gatk, coverage
---


# Question
Title: Error when trying to move multiple file types in wdl
I tried to add a "gather\_results" task to a workflow that takes in an array of Files and moves output files to a new folder of a specified name.

`Cannot interpolate Array[File?]+ into a command string with attribute set wdl PlaceholderAttributeSet(None,None,None,Some( ))]`

The wdl file passes womtools and I have no idea what a PlaceholderAttributeSet is, which makes debugging quite difficult.

--

**gather files task**

```
task gather_result{
  input{
    String dir
    String? subDir
    Array[File?]+ files
  }

  command
  {
    mkdir --parents ${dir}/${default="" subDir}; /bin/mv -t ${dir}/${default="" subDir} ${sep=" " files}
  }
  runtime
  {
    docker_image: "pipeline"
    memory: "2 GB"
  }
  output
  {
    String out = stdout()
  }
}

```

How it's called:

```
  call gather_result as gather_files {
      input: files=[create_hic.output_hic[1],
                    dedup.deduped_bam,
                    calculate_stats.stats[0],
                    calculate_stats.stats_hists[0],
                    calculate_stats.stats_json[0]],
            dir=OutputDir
    }

```

# Answer
> 1 votes
I think you'll need to convert your `Array[File?]+` to a `Array[File]` before you can use `~{sep=" " files}` in the command section.

You can use the `select_all` function for that, which will only keep valid files.

https://github.com/openwdl/wdl/blob/main/versions/development/SPEC.md#arrayx-select\_allarrayx

---
Tags: wdl
---


# Question
Title: How to download the whole genome of an organism from NCBI using R or Python?
I need to download some of the genomes of human gut microbiome in FASTA format. Anybody has any script to download?

# Answer
I'd suggest the `bit-dl-ncbi-assemblies` from the Bioinformatics Tools package.

```
$ bit-dl-ncbi-assemblies

 --------------------------------  HELP INFO  ---------------------------------

  This program downloads assembly files for NCBI genomes. It takes as input
  assembly accessions (either GCA_* or GCF_*) and optionally a specification of
  which format to download (currently: genbank, fasta, protein, or gff). For
  version info, run `bit-version`.

    Required input:

        - [-w <file>] single-column file of NCBI assembly accessions

    Optional arguments include:

        - [-f <str>] default: genbank
                  Specify the desired format. Available options currently
                  include: genbank, fasta, protein, gff, feature_tab, report, stats.

        - [-j <int> ] default: 1
                  The number of downloads you'd like to run in parallel. Write speeds
                  can become problematic generally with around 15 or more, so a max
                  of 12 will be used to be safe even if more are requested.

    Example usage:

     bit-dl-ncbi-assemblies -w ncbi_accessions.txt -f protein -j 4

```

Alternatively, you can use this code from the Biopython cookbook:

```
import os
from Bio import SeqIO
from Bio import Entrez

Entrez.email = "A.N.Other@example.com"  # Always tell NCBI who you are

seq_ids = ["EU490707"] # Your accessions go in this list
for sid in seq_ids:
    filename = f"{sid}.fasta"
    if not os.path.isfile(filename):
        # Downloading...
        net_handle = Entrez.efetch(
            db="nucleotide", id=sid, rettype="fasta", retmode="text"
        )
        out_handle = open(filename, "w")
        out_handle.write(net_handle.read())
        out_handle.close()
        net_handle.close()
        print("Saved")

```
> 3 votes

---
Tags: computational-biochemistry
---


# Question
Title: GSEA preranked different results opposing phenotypes
I am running an analysis of two phenotypes: mutant and control.

I have a statistic used to rank the different genes, and the ranking is exactly the opposite between mutant and control. Moreover, there are no duplicates in the ranking (which excludes that genes are ranked alphabetically). I uploaded the test file here.

My expectation is that running GSEA using mutant as positive phenotype, should give me the same results (but in opposite direction) than running GSEA using control as a positive phenotype - since there are only two phenotypes, and the ranking statistic is computed mutant vs control.

However, when I run Gene Set Enrichment Analysis using mutant as positive phenotype, the results I get are not the same with respect to running GSEA using mutant as negative phenotype.

I initially thought that it might be related to the choice of the random seed, however, I get the same results even setting the same random seed in all runs.

The following is an example of an enriched result of the positive phenotype:

The same results with the negative phenotype:

The ranking is exactly the same, however, the leading edge is different in the two cases. It's probably due to the fact that the ranking sum goes in one direction (increasing the enrichment score) in both phenotypes. This is probably the cause of different scores in the opposite phenotypes, however, which is the result I should use in this case?

# Answer
There are two things that I notice from these plots that are unusual for GSEA analyses that people do:

1. The zero cross is at the extreme end of the plot, meaning that *all* selected genes have a positive test statistic. I think it's a good idea to separate out the positive and negative test statistics (because then the zero point cross is at zero on the Y axis), but other people disagree with me, preferring that the plots show *all* genes, rather than positive or negative genes. In any case, GSEA seems to work best when the ranked gene list is unfiltered.
2. The "negative phenotype" plot has an unusual shape, with high enrichment score boosts for what should be low-ranked genes (and, now that I look at it, the first "positive phenotype" plot has similar issues). Most often, I see GSEA plots that have gradually decreasing enrichment score boosts as the graph moves further down the list of genes. The "negative phenotype" plot also doesn't have any sign changes; the bottom bits should be blue, not red, if the flip has been done correctly.

I suspect what's happened here is something like that there is no test statistic in the list of ranked genes, and the gene rank is being used as a distinguishing statistic. For the GSEA to look "normal", the gene list should have a field that indicates signed difference (e.g. mean expression of mutant minus mean expression of control). With a list presented in this fashion with signed difference, it doesn't really make sense to compare mutant and control phenotypes separately, as the test statistic incorporates both at once.
> 1 votes

# Answer
What do you mean by having the same results but **opposite** direction?

What GSEA does is that it goes through the gene list from the top to bottom, whenever it encounters a gene that belongs to a gene set, it adds a positive score to that gene set, and if not, it penalizes the score. That is what you see in the green line in the results.

Now imagine that you reverse the order of your input gene list, it is obvious the results are going to be different if the program still goes through the list from top to bottom.

I think as a general rule, whenever you change the order of input gene list, results are going to be different.
> 0 votes

---
Tags: gsea
---


# Question
Title: How to download data from SRA in Linux systems via the command line?
My workflow for downloading data from SRA has been the following:

1. Access SRA Run Selector.
2. Enter the accession number for the project of interest.
3. Download "Accession List" for the "Total" data. The accession list is basically, as the name says, a `.txt` file listing all accession numbers for each run in a given bioproject.
4. For each accession number in the accession list I run the following command:

```
$ fastq-dump --accession <accesion_number> --outdir .

```

How to download data from SRA in Linux systems via the command line? What are fast and scalable ways of achieving this? What should one know before downloading data?

Although there are some similar discussions in this StackExchange (eg. Trying to retrieve sra data using sra explorer, Downloading dataset from SRA (SOLiD Platform), Downloading SRA Files from AWS, etc).

The discussion on What's the best way to download data from the SRA? Is it really this slow? seems like a good starting point, but I wasn't able to find a straightforward guide on how to download data from SRA, e believe that a well thought and complete answer would help beginners alike. I am preparing an answer to the problem myself, but I thought that the community would probably give some interesting insights (:

# Answer
So far, **grabseqs** is the easiest option, a wrapper for `fastq_dump` and `fasterq-dump`, you can install it with conda, I recommend you to use an environment:

```
conda create -n grabseqs -c louiejtaylor -c bioconda -c conda-forge grabseqs

conda activate grabseqs

grabseqs sra SRRNNNNNNNN

conda deactivate

```

And download directly projects, runs, and bioprojects from SRA with `grabseqs sra SRRNNNNNNNN`. Check it out here: https://github.com/louiejtaylor/grabseqs

You can install in a different way just in case the first one doesn't work:

```
conda create -n grabseqs -yc bioconda parallel-fastq-dump pigz
conda activate grabseqs
pip install grabseqs

grabseqs sra SRRNNNNNNNN

conda deactivate

```
> 2 votes

---
Tags: genomics, public-databases, transcriptome, data-download
---


# Question
Title: programmatically select and label peptides within pdb file
I have a pdb file of a protein structure, and I would like to programmatically select certain peptides and label them in an annotated pdb file output, so that when the output pdb file is opened again by one of the 3D viewers, one can highlight the annotated peptides.

I can do this interactively using UCSF Chimera or https://www.rcsb.org/3d-view,

and with Chimera, I can even save a session which is a python file that contains the lines of my manual annotations as lists of coordinates.

whereas in RCSB, saving a session produces a molx compressed file that includes hashes such as this one, where `myseq` is defined by a range within the file:

```
              {
                "parent": "yNaD8QxiUNHaHCvxrC5rew",
                "transformer": "ms-plugin.structure-component",
                "params": {
                  "type": {
                    "name": "bundle",
                    "params": {
                      "hash": 1673297064,
                      "elements": [
                        {
                          "groupedUnits": [
                            [
                              0
                            ]
                          ],
                          "set": [],
                          "ranges": [
                            1669,
                            1777
                          ]
                        }
                      ]
                    }
                  },
                  "nullIfEmpty": true,
                  "label": "myseq2"
                },

```

So I want to do the same programmatically, e.g.:

```
annotate_sequences --input file.pdb --sequences 1:ARNDC,2:CDNRA,foo:VYWTSP

```

I expect an output file where these three peptides are annotated with the labels `1`, `2`, and `foo`, so when someone opens the file in a 3D structure viewer, then can select/deselect the annotations.

Any ideas on how to do this?

# Answer
> 2 votes
There is no consensus on how different programs do representations (example of NGL vs. PyMOL) or what they are called so there is not official way to store selections.

Annotations can be `REMARK` lines, which have a specific format or simply lines after a hash at the bottom of the file get ignored by most viewers but that specific program may read. Here is a scored Rosetta PDB file for example:

```
HEADER                                            14-JUN-21   XXXX
EXPDTA    THEORETICAL MODEL
REMARK 220
REMARK 220 EXPERIMENTAL DETAILS
REMARK 220  EXPERIMENT TYPE                : THEORETICAL MODELLING
REMARK 220  DATE OF DATA COLLECTION        : 14-JUN-21
REMARK 220
REMARK 220 REMARK: MODEL GENERATED BY ROSETTA
REMARK 220  VERSION 2020.49+release.fac69dd
ATOM      1  N   MET A   1     107.793 134.398   6.022  1.00  0.00           N
ATOM      2  CA  MET A   1     108.468 134.162   7.299  1.00  0.00           C
ATOM      3  C   MET A   1     108.644 135.405   8.219  1.00  0.00           C
ATOM      4  O   MET A   1     108.050 135.416   9.303  1.00  0.00           O
ATOM      5  CB  MET A   1     109.782 133.398   7.080  1.00  0.00           C
ATOM      6  CG  MET A   1     110.650 133.262   8.335  1.00  0.00           C
ATOM      7  SD  MET A   1     112.224 132.454   8.016  1.00  0.00           S
...
# All scores below are weighted scores, not raw scores.
#BEGIN_POSE_ENERGIES_TABLE
label fa_atr fa_rep fa_sol fa_intra_rep fa_intra_sol_xover4 lk_ball_wtd fa_elec pro_close hbond_sr_bb hbond_lr_bb hbond_bb_sc hbond_sc dslf_fa13 omega fa_dun p_aa_pp yhh_planarity ref rama_prepro total
weights 1 0.55 1 0.005 1 1 1 1.25 1 1 1 1 1.25 0.4 0.7 0.6 0.625 1 0.45 NA
pose -17624.7 2135.24 10828.2 38.8161 591.762 -521.596 -5634.4 9.63107 -1337.76 -585.258 -616.219 -453.443 0 197.59 3291.97 -602.3 2.74143 1088.19 99.1601 -9092.4
MET:NtermProteinFull_1 -2.63316 0.99093 2.70981 0.01206 0.04738 -0.17936 0.08728 0.00151 0 0 0 0 0 0.35192 1.9909 0 0 1.65735 0 5.0366

```

### Hacked b-factors

Another way is to store numerical data as b-factors. The AlphaFold2 data for example stores pLDDT in the b-factor column, while consurf stores conservation score there. So a boolean 1=show, 0=hide could be done for example, which would allow one to select these, e.g. in PyMOL:

```
# store chosen selection...
alter chosen, "b=1"
sort
save foo.pdb

# load it back...
select chosen, b<0

```

Although this makes no sense as one can save a session file... Also note, that for multiple selections, you'd need create different models, but whatever happens the other non-coordinate column, occupancy, should not be touched.

### Sharing

If the problem is that you want to share a protein with someone and want to let the recipient know about features of interest but are using different software than you, you may want to consider Michelanglo, a website I wrote that allows one to share a webpage with an interactive protein viewer (using NGL.js, hence the name) and a side description with custom annotations that control the protein representation (example).

---
Tags: protein-structure, pdb
---


# Question
Title: Converting Fastq files to Fasta files on Ubuntu
I am new to bioinformatics and programming. I tried to convert my fastq files to fasta file, but got this. What could be wrong?

```
(base) jszheng@Greater:~/seqtk$ seqtk seq -a KF0124_1_1.fastq > KF0124_1_1.fasta
[E::stk_seq] failed to open the input file/stream.

```

# Answer
*Answer from @haci, converted from comment:*

You should run your command within that directory (the directory where the fastq file is), this does not seem to be the case as far as I see. Alternatively you can provide the absolute path of your fastq file.

`cd` into the directory where the fastq file is (you should be able to see your file of interest when you type `ls`) and then run the command in your question:

```
cd ./genome_diirt
seqtk seq -a KF0124_1_1.fastq > KF0124_1_1.fasta

```

Alternatively you can provide the absolute path of your fastq file:

```
seqtk seq -a ./genome_diirt/KF0124_1_1.fastq > ./genome_diirt/KF0124_1_1.fasta

```
> 2 votes

---
Tags: fasta, genome, sequence-analysis, fastq, linux
---


# Question
Title: DESeq2. Which is better: using the altHypothesis argument in the results function or manually filtering for your desired P-value and Log2FoldChange?
In DESeq2 you can identify your differentially expressed (DE) genes using the results function. I noticed people in my lab using the results() function with the minimum number of arguments supplied (below) and then converting the resulting object to a data frame filtering for the genes they considered DE (where they considered a fold change of +/- 2 and a p-value of 0.05 to be DE).

```
res1 <- results(dds, contrast(condition, treated, untreated))

```

Alternatively, I noticed that the results function allows you to set both the P-value and minimum fold change when you use the 'altHypothesis' argument (below). The resulting list from using this method is considerably shorter than the list from manually filtering for the same P-value and fold change. I am aware that that by using the 'altHypothesis' argument it comparisons are made using a Wald test rather than a LRT test, however, I am completely unfamiliar with these two tests and I was not able to understand their significant differences in the context of RNA-seq analysis.

```
res1 <- results(dds, contrast(condition, treated, untreated), alpha = 0.1, lfcThreshold = 1, altHypothesis = "greaterAbs")

```

Due to the difference in size of the lists, and the goal of being consistent moving forward for the lab, does anyone have any insight as too which one is better? Neither myself or my lab are particularly used to this analysis but we would like to try to follow best practices as much as possible.

Thank you

# Answer
If you have an existing understanding of a fold change that you would consider to be biologically relevant, it will be better to use the lfcThreshold for comparisons. Working out that fold change is the hard bit.

When using the Wald test, the p-value is associated specifically with that threshold, rather than the difference from zero.

This is explained in the DESeq2 Analysis Vignette:

> It is also possible to provide thresholds for constructing Wald tests of significance. Two arguments to the *results* function allow for threshold-based Wald tests: `lfcThreshold`, which takes a numeric of a non-negative threshold value, and `altHypothesis`, which specifies the kind of test. Note that the alternative hypothesis is specified by the user, i.e. those genes which the user is interested in finding, and the test provides $p$ values for the null hypothesis, the complement of the set defined by the alternative. The `altHypothesis` argument can take one of the following four values, where $β$ is the log2 fold change specified by the `name` argument, and $x$ is the `lfcThreshold`.
>
> * greaterAbs - $|β|\>x$ \- tests are two-tailed
> * lessAbs - $|β|\<x$ \- p values are the maximum of the upper and lower tests
> * greater - $β\>x$
> * less - $β\<−x$
> 1 votes

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


# Question
Title: SCTransform Warning: in theta.ml(y = y, mu = fit$fitted) : iteration limit reached
Running SCTransform on my seurat object produces the warning:

```
Warning in theta.ml(y = y, mu = fit$fitted) : iteration limit reached

```

What does this warning mean? Can I trust the output?

SCTransform runs through the data producing this warning until the progress bar reaches 100%. The resulting seurat object **does** have an SCT assay and I have successfully done downstream analysis with it. I am concerned that my SCTransform has not successfully done what it is supposed to do. Is there a way to cross-check?

# Answer
The explanation can be found in link posted by **ATpoint** in his comment:

> I get at least 10 posts on GitHub, GitMemory and Biostars just on the top page.

The first author of the SCTransform paper, Christoph Hafemeister had posted this answer in a similar query on the satijalab GitHub page https://github.com/satijalab/seurat/issues/1378

"These warnings are related to the estimation of the theta parameter of the Negative Binomial regression model used to describe the molecule count distribution of each gene. When a gene has very few non-zero observations estimation of that parameter is not reliable and that's what the message indicates. This is not a problem, since the estimated theta parameters are not used directly for normalization but are regularized (sharing information across genes and thus lowering uncertainty)."
> 2 votes

---
Tags: scrnaseq, seurat, sctransform
---


# Question
Title: Where can I find a basic example for testing association between a phenotype and each of the variants of a chromosome?
I am a statistics student trying to find data suitable for a particular method. From what I have read so far, it is of interest to determine if variants within a chromosome are associated with a phenotype. I would like to understand this more and if possible, find an example of this application. In particular, what test statistic is used when associating each chromosome variant to a phenotype? Also, is it possible to gain access to examples of SNP Chip data fitting this context?

# Answer
There are a very large number of such methods. I would recommend googling "GWAS" or "genome wide association analysis" and you will find many many results.

One relatively simple example would be this plant genotype-phenotype dataset, associated with this publication. Some basic summary description of the GWAS dataset can be found here. By now the methods in this paper are very old and naive seeming, but it gives an idea for how people approached this problem in the golden age of GWAS. The data are small and simple and full of strong associations. Additionally, the phenotype measurements are based on isogenic inbred lines, which are more or less homozygous, making the analysis much easier, interpretation a lot simpler.

Another point in favor of using the plant data is that there is actually some strong genotype-phenotype association with a really strong mechanistic basis, like the various flowering time associations. For better or for worse, such unambiguous associations are often absent in e.g. human phenotypes of interest.

For some background on the general problem of controlling for population structure in GWAS, see the old EMMA paper. A much more recent example of current statistical methodology that has been making the rounds on Twitter can be found here.

For very large catalogs of GWAS data in humans, I can point you here or here. The UK Biobank dataset in particular is very commonly used and quite high quality.
> 0 votes

---
Tags: statistics, database, public-databases, snp-chip
---


title: Parallelizing microRNA targets

text:
# Question
Title: Parallelizing microRNA targets
I am trying to look for miRNA targets using a file called `Conservedfamily.txt` from the Zebrafish target scan fish website. I have written a python program to extract each of the miRNAs which works fine. However, for a big list of microRNAs, I must provide a name every time in line 22 and get one individual output file. I am thinking is there any option to parallelize the process so that somehow, I can provide a list of all microRNAs and obtain all my targets at the same time. The python code is attached below. I have posted the question in biostars but go no reply from them for 4 days. Hence I am posting the question in Stack exchange. The python code is attached below.

```
#!/usr/bin/python

#Import relevant modules

import sys

conservedfamily = open(sys.argv[1], 'r')
targetgenelist = open(sys.argv[2], 'w')
targetgeneid = open(sys.argv[3], 'w')

geneid = []

targetgenelist.write('miR family' + '\t' + 'Gene ID' + '\t' + 'Gene Symbol' + '\t' + 'Transcript ID' + '\t' + 'Species ID' + '\t' + 'UTR start' + '\t' + 'UTR end' + '\t' + 'MSA start' + '\t' + 'MSA end' + '\t' + 'Seed match' + '\n')

for line in conservedfamily:
  splitline = line.strip().split('\t')

  if splitline[0] == 'miR-210':
    targetgenelist.write(line)
    if not splitline[1] in geneid:
      geneid.append(splitline[1])
      targetgeneid.write(splitline[1] + '\n')

conservedfamily.close()
targetgenelist.close()
targetgeneid.close()

```

```
python miR99target.py 'Conserved_Family_Info.txt' 'miR-210targetlist.txt' 'miR-210targetid.txt'

```

# Answer
> 2 votes
I'd suggest loading in all of the miRNA tags as a `set()`. Then, you can just check the set membership each time.

```
#!/usr/bin/python

#Import relevant modules

import sys

conservedfamily = open(sys.argv[1], 'r')
targetgenelist = open(sys.argv[2], 'w')
targetgeneid = open(sys.argv[3], 'w')

# Create a set of all miRNA tags. Each tag should be on its own line.
with open(sys.argv[4], 'r') as miRNAlist:
    miRNA_set = {line.strip() for line in miRNAlist}

geneid = []

targetgenelist.write('miR family' + '\t' + 'Gene ID' + '\t' + 'Gene Symbol' + '\t' + 'Transcript ID' + '\t' + 'Species ID' + '\t' + 'UTR start' + '\t' + 'UTR end' + '\t' + 'MSA start' + '\t' + 'MSA end' + '\t' + 'Seed match' + '\n')

for line in conservedfamily:
  splitline = line.strip().split('\t')

  # Check if current tag is in the desired set
  if splitline[0] in miRNA_set:
    targetgenelist.write(line)
    if not splitline[1] in geneid:
      geneid.append(splitline[1])
      targetgeneid.write(splitline[1] + '\n')

conservedfamily.close()
targetgenelist.close()
targetgeneid.close()

```

One additional comment: I'd suggest using the python module for csv reading/writing, even for quick and dirty scripts as well as using context managers for reading/writing files.

---
Tags: rna-seq, genomics
---


# Question
Title: Extracting p-values using corrr package in R
I'm performing correlation analysis using the corrr package in R, but I'm unable to extract the p values. Does anyone know how to get the p values in corrr package in R?

# Answer
> 5 votes
I would suggest looking at some of the documentation.

As far as I can tell, the intention of the package is to visualize and organize correlation coefficient estimates. It does not seem that the author of the package has any interest in p-value estimation.

If you are interested in computing p-values, I would suggest using a different tool, such as the builtin `cor.test()` function.

# Answer
> 1 votes
Here is a suggested workaround:

* Generate a p-value table using one of these tools:
    + `RcmdrMisc::rcorr.adjust()`
    + `psych::corr.test()`

Then use `corrr:as_cordf()` to generate the tidy dataframe with the same structure as a regular corrr correlation table.

Personally, I have used the `Hmisc::rcorr()` package, and it works well:

```
Hmisc::rcorr() %>% `[[`('P') %>% corrr::as_cordf()

```

---
Tags: r, correlation, p-values
---


# Question
Title: Tool for rna/lna melting temperature prediction
Is there any tool available to predict the melting temperature for rna/lna oligos with its perfect complement ?

I know this amazing website from qiagen/geneglobe that estimates the tm value.

But is there any command line tool that does this?

I can estimate the normal tm value with RNAcofold and the -T option but there is no lna option.

I searched for a while but i can only find online services.

# Answer
There are several programs for simulating/predicting DNA melting but you could try them for RNA. You can look here http://bioinformatics.org/meltsim/wiki/Main/HomePage.
> 1 votes

# Answer
Both Bio.SeqUtils.MeltingTemp of Biopython and Primer3 offers different methods/algorithms to calculate the Tm value.
> 1 votes

---
Tags: rna
---


title: Help with Seurat QC ambiguity

text:
# Question
Title: Help with Seurat QC ambiguity
I have ***four PBMC samples from 10X scRNA-seq***

```
> cancer
An object of class Seurat
36601 features across 18338 samples within 1 assay
Active assay: RNA (36601 features, 0 variable features)
>

```

And this is Seurat QC plot

If I am not wrong, samples were of relatively low quality, with gene expression data revealing the presence of mitochondrial genes, as well as MALAT-1, which are suggestive of poor sample quality (dead/dying cells).

Any way by this plot I filtered cells to remove data that had more than `20 "percent" mitochondrial expression` , `> 2000 features` or `< 100 features`

```
cancer <- subset(cancer, subset = nFeature_RNA > 100 & nFeature_RNA < 2000 & percent.mt < 20)

```

By this I lost most of cells

```
> cancer
    An object of class Seurat
    36601 features across 6,883 samples within 1 assay
    Active assay: RNA (36601 features, 0 variable features)
    >

```

Please, can somebody experienced in scRNA-seq inspect this plot and my thresholds and tell me if I am wrong, any suggestion, if losing this number of cells is normal

Thank you

# Answer
> 1 votes
*Answer from @haci, converted from comment:*

I think you should not be using a cut-off higher than 20 or 25 for the reasons I explained above.

I would advise against using cells with more than 20-25% percent.mito. This metric is used as a proxy for cell death and there is not much point keeping these cells if the cell death is so severe to the point where let's say half the reads map to the mitochondrial genome (with a 50% cut-off). There will be many caveats one being apoptotic cells from different lineages clustering together, not interesting unless you are studying apoptosis related phenomena. In my opinion 6000 good quality cells are better than 10000 apoptotic cells.

# Answer
> 1 votes
From the graphs, you're losing most of the cells through the `percent.mitochondrial` threshold. The distribution is so broad on that mitochondrial plot that I probably wouldn't filter for mitochondrial expression (except possibly at ~95%). If I did have a lower threshold, I'd put it at around 50% rather than 20%, because that's where there's a change in density on the mitochondrial plot.

This suggests that there's a very large mitochondria-rich population in the PBMC dataset. I don't have a good understanding of biology and cell populations, but if I were presented with these data I'd probably try clipping out those cells with 65-95% mtDNA and see if there is a strong enrichment for markers that correspond to any particular cell population (e.g. platelets).

---
Tags: scrnaseq, seurat, 10x-genomics, qc
---


# Question
Title: How are tertiary protein structures represented in computers?
What data structure is used in representing protein structures in computers so that we can apply algorithms?

Matrix or Graph or tree?

**P.S.** I am absolutely new to Bioinformatics.

# Answer
> 4 votes
> comments from above fleshed out.

Both: A graph network is a handy way to represent a molecule, but a matrix makes for fast maths.

A graph network is a list of `nodes` (atoms), a list of `edges`, the connections between nodes and a list of `properties` associated with the nodes (symbol, name, index, residue index, residue number, chain, coordinates etc.) and occasionally a list of those of the edges (single, double etc.).

A protein is a big molecule. Same rules apply. A small compound can be represented with a string, i.e. SMILES, which is a form of graph network even if it looks nothing like one.

In a mol file and its sdf/mol2 derivatives, a tabular representation is used to store coordinates (`properties`) of the atoms of the small molecule (top part below), while the connectivity has an edge-node relationship, i.e. are edges. Note that the index of the properties in the first block is the id of the atom.

```
     RDKit          3D

 12 12  0  0  0  0  0  0  0  0999 V2000
   -5.5310   24.4560   -8.2060 O   0  0  0  0  0  0  0  0  0  0  0  0
    0.9560   21.4190   -7.6300 O   0  0  0  0  0  0  0  0  0  0  0  0
   -5.4170   23.9250   -6.8750 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.2770   19.3290   -8.7330 O   0  0  0  0  0  0  0  0  0  0  0  0
   -3.9460   23.5670   -6.5510 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.3270   19.7590   -6.4400 O   0  0  0  0  0  0  0  0  0  0  0  0
   -3.4590   22.5590   -7.5870 N   0  0  0  0  0  0  0  0  0  0  0  0
   -3.5890   23.1360   -8.9920 C   0  0  0  0  0  0  0  0  0  0  0  0
   -5.0600   23.5190   -9.2240 C   0  0  0  0  0  0  0  0  0  0  0  0
   -2.0390   22.1150   -7.2780 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.6090   21.0030   -8.2230 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.0430   20.3530   -7.7350 S   0  0  1  0  0  0  0  0  0  0  0  0
  1  3  1  0
  1  9  1  0
  3  5  1  0
  9  8  1  0
 12  2  1  1
 12  4  1  0
 12  6  1  0
 12 11  1  0
  5  7  1  0
  7  8  1  0
  7 10  1  0
 10 11  1  0
M  END

```

A PDB, CIF, etc. file is the same. It has a per-atom tabular representation of the properties, with an optional connectivity part (`CONECT` entries) for non-standard residues (not AAs and DNA).

Now when these are loaded into memory, they become an array of atoms and bonds, each with their properties —graph network.

However, the coordinates within get often converted into a matrix for speedy operations. For example to calculate an atom to atom distance matrix, converting the coordinates with the array of atoms to a matrix can be very easily be converted into a distance matrix via simple matrix calculations.

So both matrix representations and graph network representations are used.

For 3D visualisation, a pre-existing library is used and bonds and atoms become primitive shapes, whose splines may be connected to each other (in knots/vertices) or simply the bonds shapes clip through the atoms. That is there is not graph relationship.

---
Tags: protein-structure
---


title: Seurat heatmap for two conditions

text:
# Question
Title: Seurat heatmap for two conditions
I have a Seurat object of four cancers and four controls

```
> colnames(pbmc@meta.data)
 [1] "orig.ident"       "nCount_RNA"       "nFeature_RNA"     "samples"
 [5] "Condition"        "Gender"           "percent.mt"       "nCount_SCT"
 [9] "nFeature_SCT"     "SCT_snn_res.0.75" "seurat_clusters"  "SCT_snn_res.0.5"
> unique(pbmc@meta.data$Condition)
[1] "control" "cancer"
>

```

I want to have a heat map of differentially expressed between cancers and controls like `figure 3A` in this paper

https://www.frontiersin.org/articles/10.3389/fgene.2020.610682/full

Do you have any idea please how I can do that? I've looked at the integration tutorial, but can't find the heatmap in it.

# Answer
> 1 votes
See the `DoHeatmap` function in Seurat, which seems to be what that paper has used:

https://satijalab.org/seurat/articles/pbmc3k\_tutorial.html#cb50

```
> DoHeatmap(pbmc, features = myGeneList, group.by="Condition")

```

---
Tags: scrnaseq, seurat, differential-expression, 10x-genomics, heatmap
---


# Question
Title: Is it normal to have a smaller bam file after merging bam files?
I have 2 bam files that belong to the same sample and I merged them with `samtools merge`. And after merging, I realized that the merged version is a bit smaller than the sum of the other two separate files' size. I also wonder if I merged "the fastq" files of this sample first and convert it to bam would the result be the same?

# Answer
> 1 votes
If you merge a lot of BAM files you lose the overhead of the header which depending on the size of your BAM files this can be a significant difference or not.

With FASTQ files there should be less difference (as they don't have an overall header). However, the compression ratio can change depending on how you zip your files. \[EDIT: I just tested it merging two fastq.gz that were 357MB and 562MB - the result was fastq.gz that was smaller than the sum at only 826MB\]

No, the BAM header information would be lost in the conversion from BAM -\> FASTQ -\> BAM.

---
Tags: sam, samtools
---


title: PDB: Download all biological assemblies

text:
# Question
Title: PDB: Download all biological assemblies
I was wondering if it is possible to download all biological assemblies from the PDB at once (so not the standard asymmetrical files). For example, I want a file containing the coordinates of all atoms in a trimer, instead of only the coordinates of a monomer.

I saw that you can batch download with a shell script, but I'm not sure if you can also do that for biological assemblies.

# Answer
If you go to the wwPDB ftp site https://ftp.wwpdb.org/ in your browser you can find the biological unit structures (and anything else available in the PDB). The caveat of course is that a biounit structure has to have been generated for your structure.

By way of some examples, to get all biological units you could do (note the `ftp` instead of `https` in front of the URL):

```
wget ftp://ftp.wwpdb.org/pub/pdb/data/biounit/PDB/all/*

```

or to get a version organized in directories:

```
wget ftp://ftp.wwpdb.org/pub/pdb/data/biounit/PDB/divided/*

```

Note that this would take a few days to download everything, though. After the first download, you can add a `-N` to the `wget` command to just download new files.

For getting just one file:

```
wget ftp://ftp.wwpdb.org/pub/pdb/data/biounit/PDB/divided/00/100d.pdb1.gz

```

Hope this helps.
> 3 votes

---
Tags: protein-structure, pdb
---


title: Sample versus subject in scRNA-seq

text:
# Question
Title: Sample versus subject in scRNA-seq
I am reading this paper, where the authors mention

> We apply SAUCIE to the batch correcting, denoising and clustering of an 11-million cell mass cytometry dataset with 180 samples from 40 subjects in a study of the dengue flavivirus

What exactly is the difference between a “sample” and a “subject” in scRNA-Seq?

PS: here’s a possibly relevant visual from the paper, but I still couldn’t understand the difference

# Answer
> 1 votes
A subject is an individual or donor. You can take multiple samples from the same subject, e.g. from different tissues, different timepoints etc.

---
Tags: scrnaseq
---


# Question
Title: shiny dashboard to visualise RNA seq data
I am trying to build a shiny dashboard for the visualization of different RNAseq data sets.

However, I am encountering problem with reactive input, the code is always using the second dataset that I have loaded/read in. Anyone has faced a similar issue? If yes, how did you resolve it.

My code is

```
    library(shinydashboard)
library(uwot)
library(DESeq2)
library(gridExtra)
library(tidyverse)
library(RColorBrewer)
library(DESeq2)
library(pheatmap)
library(DEGreport)
library(vsn)
library(RColorBrewer)
library("genefilter")
library(org.Hs.eg.db)
library(dplyr)
library(tidyverse)
library(fgsea)
library(clusterProfiler)
library(ggplot2)

set_1<-load("C:/Users/abn/Documents/Shiny/DashBoardTutorial/TeData2.RData")
set_2<-load("C:/Users/abn/Documents/Shiny/DashBoardTutorial/TeData1.RData")
data_list = list(set_1=set_1,set_2=set_2)

ui <- dashboardPage(
    dashboardHeader(title = "Data Visualizer", titleWidth = 300),
    dashboardSidebar(
        width = 300,
        sidebarMenu(

            menuItem("Datasets",  icon = icon("cog"),
                     selectInput("Datasets", "Datasets:", choices = list("sample1" = "set_1", "sample2" = "set_2"),
                                 selected = "set_1")),
            menuItem("Quality Control", tabName = "widgets", icon = icon("th")),
            menuItem("Differential Genes", tabName = "widgets2", icon = icon("th")),
            menuItem("Downstream", tabName = "widgets3", icon = icon("th"))

         )
    ),
    ## Body content
    dashboardBody(
        tabItems(
            # First tab content
            # Second tab content
            tabItem(tabName="widgets",
                   h2("widgets"),
                   plotOutput("widgets"),

            ),
            tabItem(tabName = "widgets2",
                    h2("Widgets2 tab content"),

            ),
            tabItem(tabName = "widgets3",
                    h2("Widgets3 tab content"),
                    plotOutput("widgets3"),

            )

        )
    )
)

server <- function(input, output) {

    datasetInput <- reactive({
        df <- data_list[[input$Datasets]]
        })

        output$widgets <- renderPlot({
            datasetInput()
            par(mfrow=c(1,2))

            boxplot(counts(dds, normalized=F), outline=F, col=dds$condition, medcol = "white", cex.axis=0.6, main="Without Normalization")
            boxplot(counts(dds, normalized=T), outline=F, col=dds$condition, medcol = "white", cex.axis=0.6, main="Normalized")

            })

        }

shinyApp(ui, server)

```

# Answer
There are two issues that I can see.

The first is that the results are loaded via `RData` files, which include information about the variable names. Assuming these work, I expect they create a variable within the R environment called `dds` which has the DESeq2 results. The second load overwrites that variable:

```
set_1<-load("C:/Users/abn/Documents/Shiny/DashBoardTutorial/TeData2.RData")
set_2<-load("C:/Users/abn/Documents/Shiny/DashBoardTutorial/TeData1.RData")

```

What should probably be done is storing the results as an R Dataset that doesn't include variable name information:

```
# create set_1 DESeq2 results called 'dds'
saveRDS(dds, "C:/Users/abn/Documents/Shiny/DashBoardTutorial/TeData2.rds")
# create set_2 DESeq2 results called 'dds'
saveRDS(dds, "C:/Users/abn/Documents/Shiny/DashBoardTutorial/TeData1.rds")

```

Then those can be loaded into the R environment using the assignment you expect:

```
set_1<-readRDS("C:/Users/abn/Documents/Shiny/DashBoardTutorial/TeData2.RData")
set_2<-readRDS("C:/Users/abn/Documents/Shiny/DashBoardTutorial/TeData1.RData")
data_list = list(set_1=set_1,set_2=set_2)

```

The second issue is that the plot output displays `dds`, rather than the `datasetInput` variable, and there is no assignment to change the definition of `dds`:

```
    output$widgets <- renderPlot({
        dds <- datasetInput() ## Note: added assignment
        par(mfrow=c(1,2))

        boxplot(counts(dds, normalized=F), outline=F, col=dds$condition, medcol = "white", cex.axis=0.6, main="Without Normalization")
        boxplot(counts(dds, normalized=T), outline=F, col=dds$condition, medcol = "white", cex.axis=0.6, main="Normalized")

        })

    }

```
> 1 votes

---
Tags: r, rna-seq
---


# Question
Title: How to use multiple .pdb files as templates for homology modelling?
I just ran an AlphaFold job for PduK from the Pdu BMC operon and got very strange results including a large circular loop of random coil. According to Mayer *et al.*, PduK should be semi-triangular protein that forms hexagonal homohexamers. To fix this issue, I am interested in performing homology modeling of PduK with other Pdu proteins (namely PduA, PduB, PduJ, PduK, PduU, and PduT, all of which have crystal structures in the Protein Data Bank) as templates. I was only able to figure out how to use one template .pdb file in I-TASSER with its "-restraint4" parameter and was not able to get MODELLER working for this purpose.

My PduK residue sequence is: VKQSLGLLEVSGLALAISCADVMAKAASITLVGLEKTNGSGWMVIKIIGDVASVQAAISTGVSFADQRDGLVAHKVISRPGDGILSHSVTPESESEPAPAPTPVVPHEEIPEDHAAPEAPQDAELISCNLCLDPACPRQKGEPRSLCLHSGKRGEA

# Answer
> 1 votes
## AlphaFold2

You are asking about **how to hybridise different threaded models**, but I would not give up on AF just yet as it can do complexes.

AlphaFold2 is conceptually amazing, but does have some issues to such an extent that I blogged about common problems I was getting so many questions. Specifically, it models one protein chain in isolation.

In the RoseTTAFold paper David Baker and co. mimicked the AlphaFold2 algorithm hinted at from the CASP14 competition. They discovered that if multiple chains were given as a protein concatenation, with suitably long polyglycine linkers, a good oligomer or hetero complex could be resolved. Quickly after the public release of the Colabs notebooks of AlphaFold2 in the Twittersphere there were reports of homo/hetero-oligomers being resolved nicely. Of particular interest is the report that if you put a spurious pairing you get bad results —i.e. a control of sorts.

This strategy many be applicable to your homohexamer.

Your long random loop could be binding something else, so if you know what you may also then want to try with that something else!

## Threading, ab initio and hybridisation

Given a template PDB and an aligned sequence a threaded model of the sequence can be made by replacing the residues where required. SwissModel, ModBase and Rosetta partial\_thread app does this in the open, while Phyre2 and I-Tasser do it hidden in the pipeline. If a span is absent, SwissModel and others will try to predict from scratch (ab initio) the span. Features of different threaded models can be stuck together (hybridised) manually or automatically as do Phyre2 and I-Tasser. The RosettaCM (comparative modelling) toolkit is really handy in that one can customise how different parts are treated by adding Rosetta constraints (restraints in crystallography parlance), ligands etc etc. It should be noted that threaded models are faithful to the template, so can be converted into the same complexes as the template, but as soon as they get hybridised they loose this.

You currently have various threaded models against a panel of homologues of your choosing and are looking to hybridise them based on your requirements —I assume to circumvent a mistake the standard I-Tasser server makes, such as messing up oligomers. So the hybridisation from RosettaCM (skipping the threading) is probably a good choice, however, requires familiarity with Rosetta.

The model would need to be "symmetric" (oligomer) before starting the hybridisation (Step 2 in the manual). But that can simply be a case of getting the oligomeric template and superimposing that many copies of your threaded model onto it. Then as your threaded model is from a different framework, it needs to be adapted to the forcefield in Rosetta which differs from I-Tasser, so needs to be minimised, using Relax.

Another minor note is that a weakness is the pairwise alignment for all of these, which for distant variants can be problematically erroneous —hence why Phyre2 and I-Tasser run a MSA and from there extract the pairwise, which is a lot more accurate. I do not know how the restraint4 parameter works.

## Edit: SwissModel

I realised just now that maybe SwissModel was not tried —it is an old threader, but does oligomers nicely. In fact, I tried your sequence and the hexamer looks decent and the internal metric is pretty good...

---
Tags: protein-structure, modelling, sequence-homology, computation
---


title: Protein fold pathfinding?

text:
# Question
Title: Protein fold pathfinding?
This is probably a long-shot since it would be so dependent upon the underlying folding algorithms. Are there known algorithms for mapping the path followed in a protein folding simulation when you have the starting conformation and the solution conformation such that you can count the number of steps involved in going from one to the other?

The starting conformation would be something generated from MSM Builder while the solution conformation would be the result of hours, days, weeks, months, etc of folding.

# Answer
> 1 votes
One way to get one of the many statistically possible paths is to start from the end and work backwards. This works with the maze puzzles on kids' menus, but also with proteins. To find how a ligand enters an active site, dynamic undocking uses tethered MD to pull a bound form out. Likewise for a protein folding starting from a folded state and lowering the attraction term of Lenard-Jones or gradually increasing the reference temperature used by the Berendsen thermostat in an MD run (as done for DNA annealing simulations).

This assumes that the system is reversible (not true for protein with a maturation step like GFP), but it is avoids getting stuck in local minima or dead-end paths in the case of maze puzzles. It also assumes that peptides start folding only once they are finished being synthesised, which is not true and actually circular permutations are an interesting case as many the permutation windows/spokes are unfolded.

However, the path down a folding funnel is more of a hike/roll down a canyon than a maze puzzles on paper as there are many possible routes to the solution. So having a bunch of trajectories makes a nice video, but does not say much...

A seemingly similar concept, but very different, is the pathway reconstruction for the evolution trajectory that leads to a multi-mutation variant of interest. This is well studied because the intermediates are not transient but have to persevere in the population until another mutation appears. If they have a too high fitness cost they will not be around long enough for a compensation mutation (cf. epistasis).

# Answer
> 0 votes
*Answer from @Llrs, converted from comment:*

I would say that the problem is defining each steps on your folding path, it is a continum, I don't know if you have seen one of those videos/gifs of a protein folding, they are quite illustrative on this aspect. Usually there isn't a clear step with a pause unless it reaches a local minimum, it could be possible to save that minimums with its folding but usually they take time and resources to calculates the folding process for a long time, so probably one ends up with just a few steps of the folding.

---
Tags: protein-structure
---


# Question
Title: What is the correct way of dealing with the analysis of data from flow cytometry?
I would like to detect a change in expression of a molecule present on a cell type by flow cytometry.

Assuming I am able to detect, using an antibody, a signal that represents the amount of the molecule I'm interested into. Also assume that the same experiment will have some noise (ideally near to zero) represented by an isotype antibody. What if the noise is somehow higher than a minimal threshold (i.e. = 1%)? How can I deal with the noise? How can I compare different samples if the noise is higher than the threshold?

I was thinking about removing those samples that have a noise above the 1% but is this ok, from a statistical point of view?

Another possibility is: can I subtract the noise from the signal and use this information for the final statistical analysis? Is this correct? Is it possible to subtract a signal from another? How is this done?

EDIT just to clarify a bit:

I have a pool of neutrophils which, in resting condition, do not express the molecule on their surface. When stimulated, these cells are supposed to show an increased amount of this molecule (again, on their surface). We would like to quantify the change in the amount of molecule (we used FACS).

We are able to quantify the quantity of surface-molecule but we also get some aspecific signal (the noise due to the isotype antiboby). What we try to do is to keep the isotype signal lower than 1% using a gate-based methodology. What if the isotype signal we get is higher than the threshold (i.e. \>1%)? The antibody sometimes work some others don't. What should we do when it returns too much noise? Should we discard the data or should we perform a normalisation over the surface-molecule using the noise?

# Answer
If your goal is to simply "to quantify the change in the amount of molecule," then you can do that at the population-average level by simply taking the average antibody signal for the different cell populations (stimulated vs. not stimulated). In that case, you don't have to throw out samples with noise above a threshold - since they are still part of the population and so should be used when calculating average.

If you want to quantify "the change in the amount of molecule" at the single-cell level, then things are a bit more complicated, and the best answer depends in part on things like the shapes of the distribution of signal from the stimulated and non-stimulated populations. And related to that, when you say "noise above the 1%," do you mean that the Coefficient of Variation (CV) is greater than 1%? Or something else?
> 1 votes

---
Tags: statistics, single-cell, flow-cytometry
---