# Question
Title: How to compact variant data to their genes without bias?
I have a dataset of genes I am trying to collect data on from public databases, to use as features in machine learning. I am trying to take some features from UCSC genome browser (e.g. number of CpG islands per gene, number of DNase clusters per gene, regulatory enrichment scores etc.) however I am not sure how to control for bias where a gene that is larger in length - and so will then have more CpG islands or higher regulatory enrichment scores simply due to gene length.

Is there a way to correct for gene length when taking/condensing variant data to individual genes?

Edit:

For reference, my machine learning model aims to predict whether a gene is the most likely to be causal for a disease (out of all the genes given to the model). The model will score the genes as a regression classification between 0 to 1 (0 being least likely to cause disease and 1 being most likely to cause disease). I plan to later further investigate the genes with the highest scores.

The model uses a variety of multi-omic features (e.g. GTEx gene expression the genes have for many tissues, GWAScatalog data, gene intolerance scores, protein-protein interaction data, drug interaction data, phenotypic scores etc.). However, I am missing epigenetic data to describe my genes so I've been looking to collect based on UCSC's variant data (CpG islands, histone modifications, DNase clusters) - however this leads to my gene length problem when I am trying to reliably take data from the variant level.

I've been plotting my features and gene length, and seen that the UCSC epigenetic data does correlate with having a larger gene length if there is a higher count of regulatory sites (0.8 r2 for some), and so this is what I'm looking to correct.

# Answer
Its very easy, just let the ML sort this out for you and that is its advantage, You're thinking of GLM style calculation where you pre-screen the data with bivariate plots, where there needs to be nice Q-Q plots and low residual.

For ML simply include the gene length as one of your parameters along with CpG etc ... and the ML regression analysis SVC, lasso, ridge, random forest will figure the relationship out between gene length and CpG. You do zero, the ML does everything, hence from a statistical point of view purists object because you don't know the relationshiop the ML has deduced between the variables, but you will get regression weights for non-DNN stuff, which will give you some idea of the impact of length.

There is the issue of transformations and that can be complicated, but I'd try untransformed data first. The only disadvantage of this approach is the user will have to input the gene size when they want to check out your training algorithm.
> 2 votes

---
Tags: gene, statistics, public-databases, genomics, variants
---


title: How to understand methylation level?

text:
# Question
Title: How to understand methylation level?
For Bisulfite Sequencing data, we can calculate **methylation level** for a specific locus, $$ \mbox{methylation level} = \frac{\mbox{# methelyated reads}}{\mbox{# total reads}} $$ Here I think, for this locus, it can only have one status, methelyated or not. Or to say, **methylation level** here can only be 0 or 1. Why does it have other values in (0,1)?

# Answer
> 6 votes
You are right that a single molecule in a single position is either methylated or not methylated. However:

First, assuming your organism of interest is diploid (or of higher ploidy) one of the chromosomes could be methylated, the other not. That would give you a level of 0.5 and can be found in imprinted regions (where the paternally inherited chromosome is methylated and maternal not, or vice versa).

Second, except if you use a pure cell population of a single cell type it is likely that there will be cellular heterogeneity in your sample, e.g. from different types of blood cells or neurons vs glial cells in the brain. Some of these may have methylated CpGs in this locus (depending on cell-type-specific gene activity) while others can be unmethylated. Depending on the relative proportions of these cell types the methylation level can be anything between 0 and 1 for a mixture of cell types.

Thirdly, sometimes CpGs get grouped together. For example, 5 nearby CpGs are combined to give a single methylation level for this locus, and then the level can be 0; 0.2; 0.4; 0.6; 0.8 and 1.

Lastly, there could be technical artifacts or technical reasons why a position was not called as methylated (insufficient quality for a single read?), leading to deviation from 0 or 1.

---
Tags: sequence-alignment, sequence-analysis, methylation
---


title: Power calculation for GWAS/EWAS

text:
# Question
Title: Power calculation for GWAS/EWAS
I want to investigate, how much sample size i needed to obtain 80% power for GWAS/EWAS studies. Phynotype trait is discrete (not case/control) for human disease.

I wonder, does anyone has came across this? Really appreciate any suggestion or any tool/software which can answer this?

# Answer
> 1 votes
Interesting question I would recommend the effective population size $N\_{e}$, and error bars around this. Ne can be calculated via $F\_{st}$ and is perhaps the most important parameter to determine the genetic population size.

Under the Fisher definition $N\_{e}$ is the number of breeding individuals in the population.

The central issue with the approach you are suggesting is that population size does not equate to genetic population size. Thus the population could be huge, but if it is genetically homogeneous it will appear very small. In contrast the number of individuals could be very small, but the genetic variation between between them very large so $N\_{e}$ would be huge.

---
Tags: bioconductor, gwas, human-genome, phylogenetics, molecular-genetics
---


# Question
Title: why in RNA seq don't we only use reference transcriptome?
I would like to ask why in RNA seq analysis (alignment step) we use sometimes reference genome instead of reference transcriptome? thank you!

# Answer
> 7 votes
Pseudoaligmnent to a reference transcriptome with something like Kallisto is perfectly valid, assuming you are in an organism whose reference transcriptome is very well documented.

But if you want to find novel transcripts or novel splicing junctions, you need to align to genome.

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


# Question
Title: why in 3' quantseq should we use reference genome and not transcriptome?
I am studying about RNA seq. I have found that in 3'prime rna seq (quantseq, tagseq) we should use a reference genome with good annotation (plus known UTRs). My question is why? Why could not we use reference transcriptome? Thank you!

# Answer
> 2 votes
It would be great if you can elaborate on where you "found" this information.

In the transcriptome annotation you have a defined number of transcript per genes, but there might be others that are present in your sample and not in the annotation. For example you have treated the cells with some extreme condition that changes the 3' processing machinery.

With 3' quantseq, you are left with the read fragment that is at the end of your transcript. If the ends of these transcripts are not part of the sequences in the annotated transcriptome, these will be unaligned.

These instances might be rare if your transcriptome annotation is good, but let's say if you find that a gene of interest which should be expressed is not, and you suspect it might be due to a different isoform used, you can always check it with the genome alignment (see where reads might be falling in the gene). With a transcriptome alignment, this will not be possible.

# Answer
> 0 votes
I'm not sure that's true:

https://www.nature.com/articles/s41598-019-55434-x

"QuantSeq. 3′ Sequencing combined with Salmon provides a fast, reliable approach for high throughput RNA expression analysis"

---
Tags: rna-seq
---


title: DNA seq cutting using Linux

text:
# Question
Title: DNA seq cutting using Linux
I need to use the cut command on a human insulin sequence row wise but i don't know the command.I already ran byte wise command but the gave me column wise results. I cant use the -f option since i need to use a delimiter with it and the seq does not have any delimiters. please answer asap

# Answer
Use either generic \*NIX `cut` (you do not need a delimiter if you use `-c` option) or `head` tools. Even better, use any of a number of bioinformatics tools for sequence manipulation, such as EMBOSS `extractseq`.

**EXAMPLES:**

```
# Write into file output.txt characters 1-20 from each line of file input.txt:
cut -c1-20 input.txt > output.txt

# Write the first 5 lines of input.txt:
head -n5 input.txt > output.txt

# It is easier if you have your sequence in fasta format.
# Write characters 1..20 and 100-110 from input.fasta to output.fasta:
extractseq input.fasta output.fasta -regions "1-20 100-110"

```

---

One of the easy ways to install EMBOSS suite is using `conda` (from bioconda channel).

**SEE ALSO:**

EMBOSS `extractseq` help: https://www.bioinformatics.nl/cgi-bin/emboss/help/extractseq
> 0 votes

---
Tags: linux
---


# Question
Title: Installation of PRANK MSA in WLS Ubuntu 20.04 LTS
I want to install PRANK on the Windows 10 Linux subsystem (Ubuntu 20.04 LTS), I have followed the installation instructions to no avail.

```
wget http://wasabiapp.org/download/prank/prank.source.140603.tgz
tar xvzf prank.source.140603.tgz
cd prank-msa/src
make
./prank

[ cp prank ~/bin ]
[ sudo cp prank /usr/bin ]

```

I followed those instructions but got stuck with `cd prank-msa/src` I don't know why they want me to enter a directory that doesn't exist. I also followed this video on how to install linux files but it wasn't helpful. Any suggestions would be appreciated.

# Answer
Their instructions are outdated and don't match what exists on their server. The problem you're running into is not at the `cd` step but even before that - at your first step even. Your `wget` fails because that file does not exist. Go to http://wasabiapp.org/download/ and see for yourself - there is no `/prank/` directory in there.

I tried using `curl`, and it wrote the 404 page as a `.tar.gz` file, so if you did that, the tar step would error out.
> 0 votes

---
Tags: sequence-alignment
---


# Question
Title: Find all the bases for given reference position
Hi im looking for the aligned bases in the reads for a given reference position. im using the following script from the pysam documentataion. I adjusted it to find the specified position. in this case 24793.

```
import pysam

samfile = pysam.AlignmentFile("generated_alignment_sorted.bam", "rb" )
for pileupcolumn in samfile.pileup("chr05_modified.copy0", 10, 52000000):
    if pileupcolumn.pos == 24793:
        print ("\ncoverage at base %s = %s" %
               (pileupcolumn.pos, pileupcolumn.n))
        for pileupread in pileupcolumn.pileups:
                print ('\tbase in read %s = %s' %
                      (pileupread.alignment.query_name.split(';')[0],
                       pileupread.alignment.query_sequence[pileupread.query_position]))

samfile.close()

```

This outputs:

```
coverage at base 24793 = 6
        base in read m419941/6207/CCS Read=419941 = C

```

For all the position I tried it doesnt print out the same number of bases as it says that the coverage is. shouldnt this print out 6 times: base in read m419941/6207/CCS Read=419941 = C for different reads?

Why is this?

# Answer
`pileupcolumn.n` shows the total number of the reads that cover this basepair (in this case 24793). However, `pileupcolumn.pileups` iterates over the reads that cover this base pair AND this base pair has a minimum quality in those reads. So it has an additional filter for the reads to be considered. You can modify the quality threshold by setting `min_base_quality` parameter when you call pileup.

To have a similar filter applied to both coverage and reads, you can use `piliecolumn.get_num_aligned()` to get the number of segments and `pileupcolumn.pileups` to get the segments themselves.
> 1 votes

---
Tags: sam, pysam
---


# Question
Title: Does the kinship and inbreeding coefficients depend on population frequency of an allele?
I am reading Section 5.2, Kinship and Inbreeding Coefficients, of Kenneth Lange, Mathematical and Statistical Methods for Genetic Analysis. There the kinship coefficient $\Phi\_{i,j}$ is defined for two relatives $i$ and $j$ as the probability that a gene selected randomly from $i$ and a gene selected randomly from the same autosomal locus of $j$ are identical by descent.

I would think $\Phi\_{i,j}$ should depend on the particular population frequency of the gene or allele at that locus. However, the book does not seem to indicate that or at least does not stress that at all. Is there a dependence or not? The same question applies to the inbreeding coefficient.

More importantly, I doubt this definition is mathematically rigorous.

---

Here is an example to make what puzzles me clearer. It must be that the probability I am computing is not what the kinship coefficient. But this is the definition means to me. Some expert please elucidate the correct meaning of the definition.

Consider the simplest pedigree tree of one family of two parents with brothers $A$ and $B$. Allele $a$ is detected at a particular locus. Let $b$ be the other allele of the same gene. We ask for the probability of the brother $B$ having allele $a$ conditioned on $A$ having allele $a$.

In the following, to clutter symbols, we abuse the symbols by use $a$ to denote the population frequency or the unconditional probability of finding allele $a$ at the locus for the whole population as $a$. The same goes for allele $b$. $a,b\in\[0,1\]$ and $a+b=1$. We solve this problem by listing as follows all possible parental gene configurations, the configuration probability (considering the symmetry of paternal and maternal loci), as well as the probability finding allele $a$ for either $A$ or $B$ at the locus.

\begin{align} \text{parent config}~~~~ & \text{Pr(config)} & \text{Pr($a$|config)} \\ aa|aa~~~~~~~~~ & ~~a^4 & 1 \\ aa|ab~~~~~~~~~ & 4a^3b & 1 \\ aa|bb~~~~~~~~~ & 2a^2b^2 & 1 \\ ab|ab~~~~~~~~~ & 4a^2b^2 & \frac34 \\ ab|bb~~~~~~~~~ & 4ab^3 & \frac12 \\ bb|bb~~~~~~~~~ & ~~~~b^4 & 0 \\ \end{align}

The joint conditional probability of both $A$ and $B$ have allele $a$ at the locus is Pr($a$|config)$^2$. So the required conditional probability is $$P(r):=\text{Pr}(B\text{ has } a|A\text{ has }a)=\frac{\sum\_\text{config} \text{Pr(config)Pr(a|config)}^2}{\sum\_\text{config} \text{Pr(config)Pr(a|config)}} = \frac{r^3+4r^2+2r+4r(\frac34)^2+\frac4{2^2}}{r^3+4r^2+2r+4r\frac34+\frac42}, \quad r:=\frac ab\in \[0,\infty).$$ $P(r)$ depends on $r$ which depends on the population frequency of the alleles as I have claimed.

Consider $\[r^3+4r^2,2r,3r, 2\]$ as a positive entry vector. $P$ is a convex combination of $v:=\big\[1,1,1,\frac34,\frac12\big\]$. So $P$ is between the minimum and maximum of the entries of $v$, and $$\frac12\le P\le 1.$$ We see the extrema are achieved at $$P(r)=\begin{cases} \frac12, \quad r=0 \\\\[0.7ex\] 1,\quad r\to\infty \end{cases}$$

# Answer
The definition as set by Kenneth Lange's book is indeed quite vague and thus unrigorous. It is a probability not conditioned on observing any particular allele at a locus thus not on the population frequency of that allele, but conditioned on the following where we focus on one particular locus. It is rigorously defined as follows

A graph of lineage contains nodes, which we call persons, of a particular locus. Each node, or person, contains exactly two genes. We draw one and only one directed edge (parental relationship) from a gene of a parent person to the gene of a child person if the latter is inherited from the former. This is a directed acyclic graph (DAG). Each gene can be connected by no more than $1$ edge directed to (as opposed to away from) it. The two genes in one person cannot be connected by edges. Two persons are said to be connected if there is an edge pointing from a gene of one node to a gene of the other node. Two persons cannot be connected by more than one edge. For a pair of connected persons and their directed edge between them, the edge is equally probable to connect to one gene as the other gene of a person, i.e. the probability is $\frac12$. A pair of genes (persons) are said to be connected if there exists an undirected (simple) path from one gene (person) to the other. Two thus connected genes are called identical by descent.

In this setting, the DAG with respect to the persons is given and deterministic whilst any directed edge in the DAG assigned to two persons are random with independent probability $\frac12$ with respect to the genes within each person.

For a pair of persons on a lineage DAG that is connected, the probability of one gene equally probably chosen amongst the two genes of one person being connected, or identical by descent, to one gene equally probably chosen amongst the two genes of the other person is the coefficient of kinship. It can be proven that the probability is $\displaystyle\sum\_{p\in P} \frac1{2^{n(p)}}$ where $P$ is the set of all distinct (simple) paths connecting one person (**not** gene) to to the other and $n(p)$ is the number of persons on person connecting path $p$.

As an additional assumption, the genes that are connected have to be the same and thus assume the same allele.

---

With regard to the conditional probability example in the question, the coefficient of kinship is not adequate to be used to compute the desired probability. Part of the reason is the question inquires about the existence of an allele rather than picking one gene with 1/2 of the probability. We need some more fundamental probability to compute that. We will demonstrate that as follows.

1. We will compute the probability of the two siblings having both allele $a$ as follows.

For a pair of full siblings, they can have exactly $0, 1, 2$ identical by descent genes. Let the respective probabilities be $r\_0, r\_1, r\_2$. Notice we are not picking genes of a person at random but talking about the existence of particular types of genes in a person. For two full siblings, some thoughts lead to $r\_0=\frac14, r\_1=\frac12, r\_2=\frac14$.

By symmetry and without the loss of generality, we fix the connecting genes in the two siblings of the respective edges. We only need to consider the connections of the edges relative to the genes of the parents.

(1) Conditioned on there is no identical by descent genes, the sibling genes can independently assume any allele. The probability of having the configuration of $ab|ab$ (not considering gene ordering) is $2^2a^2b^2$ considering the distinct ordering of the alleles, that of $aa|ab$ is $\binom41 a^3b$, and that of $aa|aa$ is $a^4$. The total probability is thus $$Pr\_0:=\frac14(4a^2b^2+4a^3b+a^4).$$

(2) When there is exactly $1$ pair of identical by descent genes, the two identical by descent genes can only assume the same allele, while the other genes can assume alleles independently. When the identical by descent genes assume allele $a$ the probability of the two siblings both having allele $a$ is $a$. When the two identical by descent genes assume allele $b$ the desired probability is then $ba^2$. The total probability is then $$Pr\_1:=\frac12(a+ba^2).$$

(3) When there are exactly $2$ pairs of identical by descent genes, each pair can only assume the same allele but independently of each other. The only way that the two siblings not to both have allele $a$ is for the two identical by descent pairs to assume allele $b$. The total desired probability for the two siblings to both have allele $a$ is then $$Pr\_2:=\frac14(1-b^2)$$

2. Now we compute the probability of any one of the siblings having allele $a$. Each gene can of course assume alleles independently. So the probability is simply $$Pr\_3:=1-b^2.$$

Therefor the desired probability of both siblings having allele $a$ given that one of them has allele $a$ is $$Pr(a):=\frac{Pr\_0+Pr\_1+Pr\_2}{Pr\_3}=\frac{\frac14a(a^3-6a^2+5a+4)}{a(2-a)}$$ since $b=1-a$. Some tedious algebra which can be carried out by, say, Mathematica shows this is identical to the expression in the question in terms of $r=\frac{a}{1-a}$. $$Pr'(a)=\frac12\Big(2-\frac1{(2-a)^2}-a\Big).$$ It is easy to see $2Pr'(a)-2$ strictly decreases for $a\in\[0,2)$ as the last two terms decreases in that interval. $Pr'(1)=0$ and thus $Pr'(a)\>0,\ \forall a\in\[0,1)$. $Pr(a)$ thus strictly increases for $a\in\[0,1\]$ and $\frac12=Pr(0)\le Pr(a)\le Pr(1)=1$.
> 2 votes

# Answer
Kinship selection has a long history in evolutionary genetics, it is where the whole "selfish gene" hypothesis was formed (Maynard Smith, Bill Hamilton, and George Price ). I would be amazed if it was not rigorous. The problem the inuitive understanding of how the co-efficient works in real data appears lacking. So for example $F\_{st}$ simulation studies demonstrate that strong values occur above 0.2.

In terms of population size impact on the calculation, this is complex because it is a consequence of $N\_{e}$, studies have been done on this - and in general extremely small populations don't adhere and tend to exhibit 'Muller's racket" mechanisms. The general opinion, e.g. in HW equilibrium, is that deviation from expected behaviour can occur due to a variety of factors not represented in the original equation, but given a reasonable $N\_{e}$, $F\_{st}$ and $F\_{is}$ the essential.

In every day terms the parameters can be explored by sensitivity analysis, so e.g. drive the genetic variation to extremes and look at the change in the coefficient, reduce the genetic variation and assess again. Look for past simulation studies etc ...
> 0 votes

---
Tags: phylogenetics
---


# Question
Title: Within and between sample count normalization
recently I came across a situation where RNASeq sample quality was expresses trough DESeq sizeFactor. So authors reported quality of their samples with respect to some publicly available datasets as good/bad if the scaling factor of their samples was equal/higher/lower than 1. What is more, they computed it using FPKM's from Cufflinks and executed geometric mean calculation over genes (as the one done in DESeq after log transforming the raw counts). Underlying logic being that with FPKM the samples were normalized "within" and with the additional DESeq, normalization "between", thus ready ready for any downstream analysis (differential expression, AS, etc.)

I haven't seen something like that before nor the use of DESeq scaling factor for estimating the quality of a sample (all samples are from the same species/tissue/condition/etc. )

Can anyone comment on why is this ok/not ok because I cannot wrap my head around the logic (then again I am just trying to replicate the study, no value contribution on my end.)

Thank you

# Answer
> 1 votes
The approach you are describing seems very strange.
Crucially, the Vignette for DESeq2 states that the model only works correctly with unnormalized counts as input:

> It is important to provide count matrices as input for DESeq2’s statistical model (Love, Huber, and Anders 2014) to hold, as only the count values allow assessing the measurement precision correctly. The DESeq2 model internally corrects for library size, so transformed or normalized values such as counts scaled by library size should not be used as input.

Normalized counts like FPKM are not supported and should be converted to counts using the appropriate import funcion.

I don't really understand the use of the normalization factors for quality control because the number of reads in a sample does not correlate well with the quality. Maybe you could link to the paper in case there is something hidden in the experimental design that would explain this.

If you are interested in quality control a useful metric in my experience is the 3' sequencing bias that can be calculated using, for example, RNA-SeQC.

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


title: Using variables with fasterq-dump?

text:
# Question
Title: Using variables with fasterq-dump?
I am trying to download multiple fastq files from the SRA NCBI database which, conveniently enough, have their IDs in a range. Building a simple loop doesn't seem to do the trick, as bash appears to interpret my variable as an additional command? Help would be appreciated. This is what I have for now:

for i in {256..345};do
    sra_code = "DRR163""$i"
fasterq-dump "$sra_code" -p -v --outdir ./output/ -e 8

done
```
</code></pre>

# Answer
> 5 votes
You have an obvious error in your shell script that should have given you an error message:

```
sra_code = "DRR163""$i"

```

Should be:

```
sra_code="DRR163""$i"

```

You cannot have spaces around the `=` sign in a variable assignment. The line `sra_code = "DRR163""$i"` means "run the command `sra_code` with the arguments `= "DRR163""$i"`. Since there is no such command, you see an error.

Here's a working version of your command:

```
for i in {256..345}; do
    sra_code="DRR163$i"
    fasterq-dump "$sra_code" -p -v --outdir ./output/ -e 8
done

```

# Answer
> 0 votes
> which, conveniently enough, have their IDs in a range

I think it's a bit risky to rely on that - for this particular case they are not going to change but you may need to download more samples in the future.

I would suggest avoiding `sra` and `fastq(er)-dump` altogether and use ena instead.

In your case, it seems you want the runs in project PRJDB7736. From here https://www.ebi.ac.uk/ena/browser/view/PRJDB7736 you can get a tabular file (see link in *Download report ... TSV*) with the sample characteristics. Column `fastq_ftp` gives you the ftp link to the fastq files that you can pass to `wget` or `curl`.

You can also get the same report programmatically with:

```
curl -s -X GET "https://www.ebi.ac.uk/ena/portal/api/filereport?accession=PRJDB7736&fields=run_accession,sample_description,library_name,library_layout,fastq_ftp&result=read_run" > ena.tsv

```

You should edit the `fields=...` list to suite your needs. The fields available are those from:

```
curl -X GET "https://www.ebi.ac.uk/ena/portal/api/returnFields?dataPortal=ena&result=read_run"

```

I don't know what you think... I find this route much more useful and transparent than sra's.

---
Tags: bash, sratoolkit, fastq-dump
---


title: Time for running ADMIXTURE analysis

text:
# Question
Title: Time for running ADMIXTURE analysis
I am trying to run the ADMIXTURE software for the first time to analyse the structure of an in-house dataset of samples comparing to 1000 genomes project ancestral populations. I am trying to determine ancestry just from the X chromosome of those individuals.

To do that, I have downloaded the X chromosome file from 1000 genomes (ALL.chrX.phase3\_shapeit2\_mvncall\_integrated\_v1b.20130502.genotypes.vcf.gz). I have removed duplicated with bcftools and converted to plink binary with the following command in plink 1.9:

```
plink --bcf ALL.chrX.phase3_shapeit2_mvncall_integrated_v1b.20130502.genotypes.bcf --keep-allele-order --vcf-idspace-to _ --const-fid --allow-extra-chr 0 --split-x b37 no-fail --make-bed

```

Next, I modified the .bim files to change the polymorphisms ID so them where the same as their position, as a way to standarize their names. Then I filtered related individuals, also with plink and found the intersection of .bim files from 1000 genomes and my in-house dataset (a bed file with microarray genotyped variants only from X non PAR region), that was already filtered to contain only variants genotyped against the plus strand (without indels).

Then I merged both datasets with the following plink command:

```
plink --bfile filename1 --bmerge filename2.bed filename2.bim filename2.fam --make-bed --out filemerged

```

That command generated a file with all the SNPS that would have more than 3 alleles in the merged file so I removed all them from both datasets and merged them again. On that merged file I pruned SNPs by linkage disequilibrium:

```
plink --bfile cromXmerged2 --indep-pairwise 50 10 0.1

```

, assign sample sex on the .fam file and filtered individual with missingnes higher than 10%:

```
plink --file mergedfile --mind 0.1 --make-bed --out

```

I also set all heterozigous calls as missing with --set-hh-missing plink command.

With the final input I runned a supervised ADMIXTURE as following:

```
./admixture finalinput.bed 3 --haploid="male:23" --supervised -j4

```

Before running of course I created the finalinput.pop file in open office calc, by creating a list with the assigned population for each sample, with my in-house dataset left blanc. I just used african, european and east asian 1000 genomes samples because my own dataset has brazilian admixed samples. The final dataset has approximately 1500 SNPs and about 1600 individuals. My computer set up is CPU: i3-2370M with two cores, 4 threads@2.4 GHz, RAM 8 GB@1333 MHz and an SSD 120 GB Sata 3.0. OS Ubuntu 20.04 64 Bits.

The thing is that ADMIXTURE have been running for three days without showing any error or producing result yet. I would like to know if it is normal to this analysis to take so long or if it could be any problem with my dataset. I would appreciate any suggestion or advice!!! Thanks in advance, Iriel

# Answer
> 2 votes
Apparently there is a bug in ADMIXTURE that prevents for converging when using at the same time the haploid mode and multithreading (flags --haploid and -j respectively). The problem is solved by avoiding using one of them.

---
Tags: microarray, bed, phylogenetics
---


# Question
Title: Need a alternative or more complex version of venn diagram in python for matching dna sequences
I am in google colab and I have combined and set up a data frame list of sequences that goes like this

Location ID Sequence

1.1 ........ A ........ AAGAGATA

1.2 ........ A........... CGGATGAT

1.3 ........ B................ ACGTGATC

1.1 ........ B ........ AAGAGATA

1.2 ........ B........... CGGATGAT

1.3 ........ C................ ACGTGATC

So I need to make a alternative to Venn diagram because there are about 20+ ID types within a couple hundred rows of list just like this. And some of the sequences repeat or are unique to its id so I need something to be like

CGGATGAT Shows up 4 times in B,C,A, and D ACGTGATC shows up 1 time in B AGTGTCTC shows up 5 times in E,A,B,C, and D.

And then have some kind of graphical visualation to show all of these common sequences and where they are located and which are unique. Using python. This might be asking for too much but I'm not sure. I would use something like vennr to make a bigger complex version of this but I don't really know R.

# Answer
To get around limitations in using Venn diagrams for set overlaps, when there are more than three sets, back around 2013 or so I created something I called an Eulergrid plot (example at the bottom of the page), which an UpsetR plot appears to recapitulate, today.

The Eulergrid code I wrote was in a mix of Perl and R; the UpsetR plot code uses R. There appears to be a Python port of UpsetR that might work, as well.
> 5 votes

---
Tags: python, ngs, phylogenetics, visualization, pcr
---


title: Unable to use phangorn::phyDat

text:
# Question
Title: Unable to use phangorn::phyDat
I am following an example from ‘Analysis of Phylogenetics Second Edition and Evolution with R’ from Emmanuel Paradis.

He is doing:

I am doing pretty much the same:

```
x=c('garras','garras','garras', 'tejido','tejido','tejido')
y <- phyDat(matrix(x), "USER", levels = unique(x))

```

but I get:

> Error in phyDat.default(data, levels = levels, return.index = return.index, : data object must contain taxa names

Screenshot:

What am I doing wrong? how come the example from the book works...

# Answer
I cannot access the book, it might be an older version or there are some steps in between. The screen shot doesn't help.

If you look at the vignette and the error, what you need to provide is a matrix of characters with defined row names, for example:

```
set.seed(111)
x = sample(letters,100,replace=TRUE)
phyDat(matrix(x),"USER",levels=unique(x))

Error in phyDat.default(data, levels = levels, return.index = return.index,  :
  data object must contain taxa names

mat = matrix(x,10,10,dimnames=list(paste0("t",1:10),NULL))
phyDat(mat,"USER",levels=unique(x))
10 sequences with 10 character and 10 different site patterns.
The states are n t s o y e c h z q m j d o h n z y x p q g u o a i h y c g t d f y t d f e u e f d i j r z x w u m a j p j g d t q y v o x p m s o z x z a y j z r u s a y k l f g l f r i l r f m t t e f x v f g q j

```
> 1 votes

---
Tags: r, phylogenetics
---


title: STRUCTURE interprets pop ID as locus

text:
# Question
Title: STRUCTURE interprets pop ID as locus
I'm using `structure` for the first time and I fail to execute a "job". For some reason, `structure` interprets the `pop ID` column as locus.

I have following `.txt` file

The first column are `sample labels`, the second column are `population IDs` and then there are 26 columns containing the locus information (i.e. 13 markers). Thus, in total 28 columns. It is a diploïd organism so each locus consists of 2 consecutive columns (as described in `structure` manual). In total there are 64 rows (so 64 samples)

I load this `.txt` file into `structure`, specifying following parameters:

-`number of individuals` 64

-`Ploidy of data` 2

-`number of loci` 13

-`missing data value` -9

-`special format: data file stores data for individuals in single line`

-`individual ID for each individual`

-`putative population origin for each individual`

The loaded file looks like:

This looks good. My samples have a separate column, `pop ID` is a separate column for possible populations an then for each `locus` there are 2 columns. Not all loci are visible on screenshot, but there are 13.

I want to calculate `Delta K` and therefore I want to run a job. However when I run it, I get following error:

I get this error for each row/individual. It appears that `structure` interprets the `pop ID` values ("Noord","Zuid", "onbekend") as 14th locus, but I have never specified a 14th locus. I fail to understand why this happens.

When I run the program without `pop ID` column, it works.

Could anyone please advise me on how to solve this error?

# Answer
> 1 votes
Hmmmm .... you don't actually need pop ID to work for Structure, providing their are no logical errors (you get a wrong answer). The bit that must work is the *a priori* number of populations in the data set. You are then going to manually (or via a script) feed a range of *a priori* values into Structure, run each and plot the likelihood. Where the likelihood plateaus thats the number of populations contained within your data. This is the key output of Structure.

Personally, I would perform a Structure run on the test data and remove the `pop ID` column and check the answers (likelihood -lnL) are identical. Thereafter use an approach such as $F\_{ST}$ to interrogate the exact members of each population delineation.

---

I think your test will work .... The only thing I can think of is the header comprises a space, This format just needs checking against the test data.

---
Tags: phylogenetics
---


# Question
Title: Get gene sequence based on the annotation
I've got the reference genome with Python like so:

```
sequences_by_chr = {}

with gzip.open("data/Zea_mays.B73_RefGen_v4.dna.toplevel.fa.gz", "rt") as f:
    for seq in SeqIO.parse(f, "fasta"):
        sequences_by_chr[seq.id] = seq.seq

```

I've also parsed a GFF3 annotation and would like to get the gene sequence based on said annotation. Suppose the following:

```
chromosome = 1
strand = "-"
start = 45286331
end = 45290076

```

Since it's on the "minus" strand, should I do the following:

`sequences_by_chr[chromosome].complement()[start:end]`

or should I use the `reverse_complement()`?

I'm completely confused by this, would appreciate any hint!

# Answer
> 2 votes
You can parse the GFF3 to generate a bed-like file with gene ranges and use that with `bedtools getfasta` to obtain a multi-FASTA file with gene sequences. An example case is shown below:

<pre><code>## get fasta and gff3 files
wget ftp://ftp.ensembl.org//pub/current_gff3/saccharomyces_cerevisiae/Saccharomyces_cerevisiae.R64-1-1.100.gff3.gz
wget ftp://ftp.ensembl.org:/pub/current_fasta/saccharomyces_cerevisiae/dna/Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa.gz

## unzip fasta
gunzip Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa.gz

## generate genes.bed
zgrep -v '^#' Saccharomyces_cerevisiae.R64-1-1.100.gff3.gz | awk 'BEGIN{FS="\t";OFS="\t"}($3=="gene"){print $1,$4-1,$5,"name",1000,$7}' > genes.bed

## create gene fasta
bedtools getfasta -fi Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa -bed genes.bed -s -fo genes.fa
```
</code></pre>

# Answer
> 0 votes
DNA sequence is always provided on the plus strand, if the gff tells you it is on the minus strand, you need to reverse complement it, because the start of the gene in your example will be 45290076 and the end will now be 45286331.

I can give you an example below, this is a potential coding sequence from yeast, no introns and on the minus strand:

```
gunzip -c test.fa.gz
>1
TCATCCCCTTATATAAGATGGGAAACGAAGGTAAAAGATTAACGATAGCAT

```

If we read it in like you did, and if you apply the complement, this is not a coding sequence.

```
chromosome="1"
start = 1
end =51
sequences_by_chr[chromosome].complement()[(start-1):end]

Seq('AGTAGGGGAATATATTCTACCCTTTGCTTCCATTTTCTAATTGCTATCGTA', SingleLetterAlphabet())

```

Getting the sequence followed by reverse complement gives you a coding sequence, starting with ATG ending with TGA:

```
sequences_by_chr["1"][(start-1):end].reverse_complement()

Seq('ATGCTATCGTTAATCTTTTACCTTCGTTTCCCATCTTATATAAGGGGATGA', SingleLetterAlphabet())

```

# Answer
> 0 votes
You should do:

```
if sequences_by_chrom[strand] == "-":
    rev_comp( sequences_by_chr[chromosome][start:end] )
else:
    sequences_by_chr[chromosome][start:end]

```

where

```
def rev_comp(sequence):
    d = {'A':'T', 'T':'A','C':'G','G':'C'}
    return ''.join([d[i] for i in sequence][::-1])

```

---
Tags: fasta, python, reference-genome, gff3
---


# Question
Title: Produce .vcf file of ALL mutations in .bam file
I got tasked with a problem, that i thought would be quite simple to solve, but turned out to be quite tricky. Our lab is running targeted mutagenesis experiments in yeast using crispr base editors. What we have done is set up an experiment where the base editor selectively mutates a region inside of a gene of interest in a culture of yeast cells. We then extract the DNA of the whole culture and perform sequencing of our gene of interest (the size of which is around 1kb). What i would like to do is plot the mutation frequency of say, C -\> G edits, across the length of the entire gene (around 300bp), with the hope that i see a spike in mutations at the site where the base editor binds. Or at least, higher mutation rates at this site compared to the background mutation rate.

These edits are quite rare (we think), and since it's in a culture of multiple cells, they don't occur at the exact same position on the gene.

I initially thought i could do this by quality filtering the reads, aligning them to the gene sequence in order to produce a .bam file, and then piping this file through variant calling algorithms such as vcftools. However, i noticed that i get very few variants out of this, even when using extremely relaxed settings (e.g. p-value \< 1 for the variant call).

**Im therefore wondering if anyone knows how to extract from a .bam file ALL the mutations across all reads that do not match the reference sequence?**

I'm aware that this will probably lead to quite a huge .vcf file, but since the reference sequence is only 300bp, it should still be manageable? It could also be the the experiment does not work, but i would at least expect some mutations in the sequencing reads compared to the reference gene sequence.

Any help would be greatly appreciated!

# Answer
Ensure that the variant caller stringency is relaxed all the way. For example in `freebayes` use `--min-alternate-count 1 --min-alternate-fraction 0` .
> 2 votes

# Answer
You could perform an mpileup with samtools. I think this would be the most sensitive approach. Here's the man page: http://www.htslib.org/doc/samtools-mpileup.html

You might want to disable BAQ correction. If you specify a reference sequence, then any letter in the fifth column would represent a non reference base. Reference bases are coded as `.` or `,`.

The reason why vcftools doesn't work is probably because it assumes a diploid individual, which is not what you have. You have lots of samples mixed together. Freebayes allows relaxation of the diploid assumption. Here is the relevant example from the documentation https://github.com/ekg/freebayes :

> Generate frequency-based calls for all variants passing input thresholds. You'd do this in the case that you didn't know the number of samples in the pool.

> freebayes -f ref.fa -F 0.01 -C 1 --pooled-continuous aln.bam \>var.vcf

You could try a cancer/normal variant caller like mutect, but I don't think you'll have much luck.

The most straightforward approach is to use samtools mpileup and to process the output yourself. You could also use bcftools mpileup to get a VCF, and then extract the allele count from the INFO field. If you are a Python lover, then you could also use pysam to both perform the pileup and process the results.
> 2 votes

---
Tags: sam, vcf, variation, variants
---


title: Software for taxonomic assignment?

text:
# Question
Title: Software for taxonomic assignment?
I have a couple of hundred bacterial sequences of 2-30 genes of interest each, recovered from metagenomics. None of them encode rRNA. Normally I'd just BLAST the one gene I already know to be reliable phylogenetic marker, where present, but the scale of my current dataset makes it a bit daunting.

Is there any software for taxonomic assignment that you could recommend? It's fine if it requires command line use. There are so many options listed online that I didn't even know where to start... Please and thank you.

# Answer
You could take a look at kraken2 It should fit your purpose.
> 4 votes

# Answer
If you like blast, maybe you can try DIAMOND \- the authors claim it's much faster than BLAST.
> 3 votes

---
Tags: blast, metagenome, taxonomy
---


# Question
Title: How to assess if this strategy is good to have some phylogenetic information?
I have several organisms and I want to calculate for each one of these their codon composition implementing a pairwise distance matrix between these codon vectors to build a phylogenetic tree. Therefore, the matrix will have on the principal diagonal 0 values (since each organism will have a distance of zero from itself) and the other matrix members will represent the euclidean distance among one organism and another one. Thus we will have a square matrix. To assess to robustness of the various phylogenetic relationships I would bootstrap this matrix building the tree on these data. Is this a good strategy?

# Answer
> 2 votes
The method you are proposing is called a distance method, because you are constructing a distance matrix and then using a clustering algorithm to make a tree. The things to watch out for are:

* Multiple substitutions, traditionally this is take care of using a Jukes-Cantor correction. You could argue this is less of a concern because this is a 64 substition matrix model rather than just 4 or 20 substition matrix for nucleotides or amino acids.
* Use the neighbour-joining clustering algorithm. Do not use a standard hierarchical clustering algorithm such as UPGMA. I'm sure R must have the nj in something like R's ape, pretty sure scipy will have something like that
* Code for bootstrap resampling. That is a method for measuring robustness

The codon model for phylogenetic trees has been implemented in MrBayes, I suspect it will be in phyloBayes. This method will be Bayesian and you should run this is parallel with your own solution.

*Concerns* Codon models are very prone to mutation saturation (which is complicated) so you must only use it for closely related sequences. If e.g. 3rd codon distances start exceeding 70% codon models go kaput, because it will select a different codon, which is a completely separate part of the matrix - completely by chance. So whilst mutation saturation in a nucleotide model (4 substitution matrix) would invalidate the tree in theory, they would be sufficient 'good' signal for key phylogenetic signals to be identified, this wouldn't happen for a codon model.

I recognise the 'concerns' might sound a bit confusing at present, but just to reiterate, codon models need sequence data that is really closely related. For this data they excel and give excellent bootstrap support.

---
Tags: phylogenetics, phylogeny, matlab, matrix
---


# Question
Title: Minimal working example for phylogenetic tree construction in R using an alignment of DNA strings and a vector of names?
I am looking at how to build a `phylogenetic tree` in `R` from **aligned** genetic information (dna sequences). There are multiple tutorials on how to go about this task eg tutorial1, tutorial2 ("ape" docs), tutorial3, and there are many more. There appears to be quite a few steps with dependencies from multiple packages: eg. "ape", "phangorn", "phytools", "geiger". I would like to learn a minimal working example (mwe) for constructing a phylogenetic tree from a basic data-structure holding dna-strings.

Eg. I have a vector of nucleotide strings, and a separate vector for the 'names' of each string in the vector, how can I construct a phylogenetic tree from this information? (not a fasta file)

# Answer
> 1 votes
A key difference from many examples provided is that they assume that the data provided is in the form of sequences rather than an already aligned sequence alignment.

```
library("Biostrings")
library("seqinr")
library("ape")
library(phylogram)
library("dendextend")

#read in the data from a format with the aligned sequences (eg msf)
msf.res <- read.alignment(file = "geneticAlignment.msf", format = "msf")
#dist.alignment and not dna.dist commonly used in tutorials
msf.res.dist.alignment = dist.alignment(msf.res, matrix = "identity")
#perform neighbor joining upon the distance set
msf.res.dist.alignment.nj = nj(msf.res.dist.alignment)
#then plot the joined sequences
plot(msf.res.dist.alignment.nj, main = "Neighbor Joining")

```

---
Tags: phylogenetics, phylogeny
---


# Question
Title: Removing Batch Effect in Heatmaps after Differential Gene Expression Analysis
I'm working on a dataset in which the first replicate of each group is one batch and the second replicate is in a second batch. After checking the `PCA` plot and seeing the batch effect in PC1, I used `removeBatchEffect` function from limma to subtract the batch effect from my count data. Then, using PCA gives me a plot that doesn't seem to have any apparent batch effect left! However, it is recommended not to use batch effect correction for differential gene analysis but use the batch variable along with the group variable in constructing the `model.matrix`. So, I did that, ran `limma/voom` on the normalized counts, and extracted differentially expressed genes. However, when I'm trying to make a heatmap from the DEGs, I still see that the samples from different batches are clustered separately, instead of seeing the clustering of replicates of the same sample. So, my question is whether I should use `removeBatchEffect` on the count data from DEGs and then use the transformed dataset for heatmaps or there is another way to fix this?

# Answer
> 5 votes
It is indeed true that for the DE analysis one should include batch into the formula to avoid changing the original counts. Still, for everything else such as plotting heatmaps use of `removeBatchEffects` is perfectly fine and (at least to me) a standard and well-accepted procedure. It essentially does not matter what you use to correct for the batch effect for the counts you use downstream. Results will probably be similar. `Combat-Seq` from the sva package is a recent adaptation of ComBat specifically for RNA-seq which (from what I understand) better deals with the integer-count nature of the data. This operates on raw counts and avoids the infamous negative values that happen at times with both `limma` and `Combat`. After applying `ComBat-Seq` on your raw counts you can normalize them as usual with edgeR (or any tool you like) and then make the heatmaps. See https://github.com/zhangyuqing/ComBat-seq. The DE results still should come from the normal DE pipeline with batch as covariate as discussed above.

---
Tags: rna-seq, edger, limma, heatmap, batch-effects
---


# Question
Title: How to modify DNA evolution model to fit actual data?
I'm working on understanding the evolution of a gene in a phylogenetic tree. I know the rates of evolution of each organism (from the tree).

I am taking a random DNA sequence with my gene, evolving it along the tree with the said rates and comparing this to the actual gene in these various organisms. I am 'evolving' the gene if randomly at a probability determined by the rate at each step in the tree.

The gene from my simulation is different from what is actually in these organisms. How do I modify my model, so that my simulation gives the same results?

# Answer
> 1 votes
Mutation is not a random process because there are many mutational biases, e.g. the transition vs. transversion ratio. Ideally use a proper evolutionary model. Maybe try a tool that incorporates substitution and indel models. I'm not up-to-date on the best packages, but e.g.

# Answer
> 1 votes
I decided that this was too involved for a comment, so I'm making it an answer.

It is still quite difficult to answer the question given the absence of methods explanation; but I think that there's some misunderstanding of phylogenetic modeling going on here, and I can at least address that. I'm going to try to break it down into its components, starting with the basic statistical idea of simulation and then going back to phylogenies.

**Statistical reasoning**

Let's say that I flip 10 coins (2 outcomes: head (H) and tails (T)), and I observe the sequence `HTHTTTHHTT`. I then fit a maximum likelihood model, and find that the best model for the data has

$Pr(H) = 0.4$,

$Pr(T) = 0.6$,

and that these two outcomes are mutually exhaustive:

$Pr(H) + Pr(T) = 1.0$

I can then use a binomial model to re-simulate data based on these parameters. I will observe various sequences, for example:

`HTTTTTTHHT` `TTTTHHTHHH` `HTHHTHHTHH`...

However, based purely on the parameters of the model, I don't expect every new simulation to recapitulate the exact same `HTHTTTHHTT` sequence, because I have reduced its information to some rate parameters for the simulation.

Each of those new sequences will yield model parameters ($Pr(H), Pr(T)$) that are similar to the training data, but they will not necessarily be particularly similar to the training data in sequence.

**Phylogenetic modeling**

Phylogenetic modeling works in the same way, such that you reduce the observations at the tips of the tree into some model with some parameters. Most of the time you are estimating e.g. a rate matrix of parameters for transition between states for a set of residues (in simplest case for DNA 4 nucleotide states and for proteins, 20 amino acid states usually). If you are working with gene sequences you are probably using the 4 nucleotide states or some derivation thereof.

Simulations usually work by starting with some fixed ancestral state at the root (which can take various states), and then evolving forward from there along each branch of the tree using the rate matrix in a continuous-time Markov chain to "replay the tape of life", by assuming that the model you've fit somewhat reflects truth.

While opinions differ somewhat on what to expect from "replaying the tape of life" (see that linked review for more), basically no one expects to get the same results as the input data, or even 90% similar to the input data. That would indicate that evolution is highly deterministic, which just isn't reliably true. Of course there are situations where such determinism happens to be the case, but they are relatively rare, and I don't think that sequence evolution is considered to be very deterministic at all.

It's a little different in that the tree defines a correlation structure of the tips, which constrains the simulation results to be somewhat more similar to the training data than we expect from coin-flipping. But there's still no reason to expect particularly similar results on the tips.

But you can also model discrete characters on a tree, which can be even simpler. For example, I used the R APE package to simulate the following character states (A or B for each tip) on a simulated tree of 10 tips with the following commands:

```
tr = rtree(10)
c_tr = rTraitDisc(tr)

```

I can then build a model for this character and this tree using the following:

```
rates = ace(c_tr, tr, type="discrete")

```

and use it to simulate new sets of states:

```
# A (index 1) is estimated to be ancestral state with 99% confidence
new_sim = rTraitDisc(tr, rate=rates$rates, states=c("A", "B"), root.value=1)
plot(tr)
tiplabels(new_sim, adj = c(1,0))

```

Look! Not a single observation of state B! Let's try it again:

```
new_sim2 = rTraitDisc(tr, rate=rates$rates, states=c("A", "B"), root.value=1)
plot(tr)
tiplabels(new_sim2, adj = c(1,0))

```

Ok, this is more similar, but there's still another state B on the tips.

**Conclusion**

Evolution is not very deterministic. The correlation structure of the tree makes the outcomes less random than independent trials as in coin-flipping, but it is not at all expected that simulation results from phylogenetic models will match the training observations well. The best that can be said is that each simulated dataset can hopefully be *reduced to the same model parameters* as the input data, but that's a very different statement from yielding the exact same outcomes.

If what you are looking for from simulation is to yield the exact same sequence states at the tips, then you are unlikely to get that from phylogenetic simulation. You might occasionally get pretty close, but that is not at all the expectation.

---
Tags: phylogenetics, simulated-data
---


title: Gap-fill assembly with PacBio reads

text:
# Question
Title: Gap-fill assembly with PacBio reads
I would like to gap-fill (or correct) a mammalian-sized assembly for which we have BACs that have been sequenced with PacBio CLR data. I have seen there are pbm22/gcpp tools available for correcting assemblies, but looking at how they've been used, they seem to correct existing assemblies built with FALCON or similar.

Can pbmm2 and gcpp be used with a reference constructed without pacbio data? E.g. given a mammalian-sized reference and a .bam file of PacBio reads for a BAC, what's the best way to gap-fill the current reference? Thx

# Answer
> 1 votes
My understanding is that PBJelly is the tool intended to solve this problem.

This person seems to have used it starting with Illumina scaffolds. I have not personally used it but have heard that it works pretty well.

Edit: note that the PacBio ecosystem is not super easy to use from scratch. However, they have made pretty good efforts to facilitate setup using conda, as for instance documented here. This person shows a conda recipe for pbjelly.

---
Tags: assembly, reference-genome, pacbio
---


title: Delete chains from a protein

text:
# Question
Title: Delete chains from a protein
We have a docking target protein with multiple chains (A,B and C). However, I just need chain A. How can I delete the 2 other chains from the docking target protein.

# Answer
> 2 votes
If the tag `r` is incorrect, the easiest way is to download and install `PyMOL`: https://pymol.org/2/

Drap and drop your pdb file and type `remove not chain A`. Then save molecule via the menu as pdb or type `save some_file_name.pdb` (will be in your home root folder).

# Answer
> 1 votes
You can also do this with pdb-tools:

```
pdb_selchain -A prot.pdb > prot_A.pdb

```

---
Tags: protein-structure
---


# Question
Title: How can I match sequences with IDs using python dictionaries?
```
0.980634         -                    572         sp|P06672|PDC_ZYMMO                         -                     568                               1.2e-263             874.5          0.0            1        1      6.6e-267  1.4e-263   874.2      0.0       5         562    15     570    12     572    0.98   Uncharacterized                         protein                   OS=Lomentospora            prolificans        OX=41688             GN=jhhlp_005678  PE=3             SV=1

```

I've asked a similar question before. I have a better understanding of python now and would like to revisit this. Can someone critique my thinking?

I have an HMM output file that contains about 20 fields. There are two sets of IDs, e.g., P006672 and jhhlp\_005678

In a separate file are FASTA-formatted sequences linked to these IDs.

Proposed solution:

1

```
   a) read in fields containing IDs, lets called them fields 8 and 6;

   b) store field values in an array or dictionary (whatever the proper data structure is for
  this). Let's called this Fields_8_and_6

```
2. a) read proteome file containing FASTA-formatted sequences

b) Grep and store as ID:sequence key:value pairs; let's call this ID\_to\_sequence
3. a) Compare Fields\_8\_and\_6 to ID\_to\_sequence

b) Write ID1:sequence1 ID2:sequence2 as pairs into one file

4. send to ClustalO

# Answer
Biopython can index a fastafile and you can access the sequences by their ID in a dict-like manner.

All that is left for you to do is extract the IDs from your file. (more Biopython)
> 3 votes

---
Tags: bioconductor, python, proteomics
---


title: B-factor analysis in nanobodies

text:
# Question
Title: B-factor analysis in nanobodies
I am doing my research on nanobodies. Currently, I am trying to do a comparison study of the flexibility of CDR loops using B-factor values. I would like to know if there is any statistical analysis tool for this or if anyone has done anything similar.

# Answer
So B-factors, for "blur factors", are how much movement an atom has in the crystal structure. They are not reliable. B-factor putty in PyMOL is the most you should push them.

## B-factors bad

This is partly a result of flexibility, but it is highly influenced by crystal packing, which is not natural. Many enzymes, for example, are seen with an active site entrance way that is closed up in a crystal structure whereas in an MD run there is a gate loop that flops around.

Additionally, a B-factor is a single value that masks lots of forms of movement. Partial occupancy after all is at the discretion of the experimenter, so you might even have a residue in two or three low energy states that result nevertheless in B-factor disorder.

Furthermore, incorrectly fitted residues will have a higher B-factor.

Also, the meaning of B-factor are mighty cryptic, and depend on the refinement —the units are Å^2 as opposed to say Å, like a standard deviation.

## MD RMSF good

The best way to determine the movement of a residue is to do an MD run. Gromacs has very nice tutorial and is the most accessible. The metric in that case is the RMSF, which is the RMSD relative to the average. Gromacs has a utility for it called rmsf. Alternatively RMSD relative to your original structure is fine —but less of a metric of movement. MD is complained about for three reasons, none are valid here:

* Tricky to set up: this is a straightforward case
* Lengthy: You would not be after large scale conformational changes so a limited amount of steps would suffice —possibly even on a laptop on a single core for a couple of hours (cluster is always better though!). Running at different temperatures and switching the intermediate steps (replica exchange) is often done to speed things up, but would just confuse the results.
* Small molecules: A nanobody has no cofactors or small molecule ligands so you'd avoid the topology parameterisation issues. If Eukaryotic made, your nanobodies may be glycosylated — play it safe and ignore.

## B-factor extraction

If you want to extract the B-factors, you can use pymol as a python module —installed with conda and not the binary/exe installer. And use python as follows:

```
pdb_filename = ""
selector = "name CA"
import pymol2
with pymol2.PyMOL() as pymol:
     pymol.cmd.load(pdb_filename)
     bfactors = {atom.resi: atom.b for atom in pymol.cmd.get_model(selector).atom}

```

From there do whatever you fancy with the b-factors. For alternative selector strings see Pymol selection algebra. One operation is to divide the values by 8 pi^2 and take the square root of that, which is the amplitude of oscillation (U) of the atom around its equilibrium position —ie. that stated in the coordinates.
> 1 votes

---
Tags: proteins, statistics, protein-structure, computational-biochemistry
---


# Question
Title: Parsing the json file from gene ontology result into dataframe
How do I parse the json file from ontology output. My file json file

```
jdata <- read_json("analysis.json", simplifyVector = TRUE)

```

I want to convert it into a `data.frame` format. I can able to take the GO term and respective p value etc. But I would like to know the genes that are getting part of the enriched list which is under the upload\_1 for example under DNA Replication there are 33 terms but 8 are from my input. So my objective is to get the 8 gene so on and so forth.

So I would like to see a `data.frame` where I want to have GO term in the first column gene mapped into the term would be second i.e upload\_\[#\] third one fold enrichment and last one p value.

This is how i would like to see this one is the text output in the column `upload_1 (450)` instead of that number i would like to see then gene names which are present in the json file.

Some progress here is my code

```
l <- fromJSON("analysis.json")

m <- lapply(
  l[[1]]$result,
  function(x) c(x$input_list$mapped_id_list['mapped_id'], x$term['id'],x$input_list['pValue'],x$term['label'])
)
m <- do.call(rbind, m)

```

output:

```
mapped_id     id           pValue       label
[1,] Character,8   "GO:0006270" 0.007540801  "DNA replication initiation"
[2,] Character,388 "GO:0009987" 3.232843e-10 "cellular process"

```

Right now I can parse out first two row.

Things I would like to append to the code is

* I would like to do for all the terms
* The first column where i am trying to get all the genes now it is coming as character while number suggests gene present. That I wont need probably since i need only gene names which is being shown as character now.

```
dput(head(b))
    structure(list(`GO biological process complete` = c("DNA replication initiation (GO:0006270)",
    "regulation of blood coagulation (GO:0030193)", "regulation of hemostasis (GO:1900046)",
    "regulation of coagulation (GO:0050818)", "platelet degranulation (GO:0002576)",
    "cell cycle G1/S phase transition (GO:0044843)"), `upload_1 (450)` = c(8,
    13, 13, 13, 16, 15), `upload_1 (fold Enrichment)` = c(11.23,
    7.09, 7, 6.62, 5.75, 5.35), `upload_1 (P-value)` = c(0.00754,
    0.000627, 0.000716, 0.00135, 0.000347, 0.00228)), row.names = c(NA,
    -6L), class = c("tbl_df", "tbl", "data.frame"))

```

Any suggestion or help would be really appreciated.

# Answer
I don't know how to do that in R but I think this is what you want:

I used http://json2table.com for view the data structure of the JSON file, then I used `json_normalize` to normalize the JSON data into a flat table using the key `[result]` to a list of DataFrames for each index in the JSON file. Then, I concatenated the list of DataFrames in a single one with `df.concat()`

```
import json
import pandas as pd

with open('analysis.json') as f:
    d = json.load(f)

dfs = []
for x, y in enumerate(d):
    data = pd.json_normalize(d[x]['result'])
    dfs.append(data)

df = pd.concat(dfs, ignore_index=True)

df.to_csv('analysis.csv', index=False)

```
> 2 votes

---
Tags: r, gene-ontology
---


# Question
Title: How can the cell line contribution be estimated from RNASeq data?
Using a laser-capture microdissection of cells a group of cells stained with the marker of interest was sequenced. In another cohort of patients (this is all human liver tissue) the whole tissue was sequenced (RNA-seq in both cases)

Can I estimate the contribution of the cells marked in the whole liver ("weight of these cells" in the liver in words of my PI)?

My gut feeling is that it can't be done this way, it would require both single cell sequencing and whole tissue sequencing to estimate the contribution of each cell line. But perhaps there is a tool that given the cell lines or the main expression of other cells lines it can be compared to using GSVA or some similar tool.

# Answer
There are couple of computational methods which try to do this (I never used them, so no experience):

1. CellMix, based on sets of marker gene lists
2. Subset Prediction from Enrichment Correlation, which is based on correlations with subset-specific genes across a set of samples.
3. Cell type enrichment, which uses our highly expressed, cell specific gene database
4. Cell type-specific significance analysis using differential gene expression for each cell type

You might have to get some reference expression levels from public databases or papers for some of the methods.

One thing to keep in mind: you cannot really compute the cells proportion, only RNA proportion. If you have a good reason to assume that RNA quantity per cell is very similar, this is a good proxy for the cells proportion in a tissue.
> 6 votes

# Answer
Cell deconvolution is mentioned in this Biostars post, which mentions CIBERSORT for immune cell mixes, and the Bioconductor package DeconRNASeq.

As far as I'm aware, it is only possible at best to get proportional representation for transcript expression from standard high-throughput sequencing results, because the sequencers and sample preparation workflow are designed in such a way that the same number of reads are output regardless of the input amount.
> 4 votes

# Answer
You can download public single-cell RNA-seq data or bulk RNA-seq for purified primary cells for reference and then use CIBERSORT or MuSiC to deconvolve. One thing to notice: to make bulk RNA-seq and single-cell reference compatible, TPM from RNA-seq can be compared with CPM in 3' (10X,Drop-seq etc.) single-cell data and TPM in full-length single-cell data (SMARTseq etc.)
> 0 votes

# Answer
Actually a month ago, I was searching for the same question and I found a few packages that can be used. First of all the functionality of them can be cell proportion estimation, cell specific gene expression estimation, or both of these functions. Secondly, their primary input, as expected, is bulk RNA-seq data and some also want a single cell RNA-seq data as secondary input. The most popular one is the child of famous cibersort, CIBERSORTX. I went through their website and read the tutorials. They have several example datasets that can be used for learning their platform. However, the two biggest flaws were that they haven't shared their source code which means we can't reproduce the result on our own systems and we have to upload all of our data to Stanford servers (security concerns!). The second problem is that there is no API access to their server and everything should be done online which means very low compatibility with other tools.

The idea, using deconvolution and deep learning to infer subpopulation of cells from bulk RNA-seq, is great. So, I searched to see if others have also done something similar and maybe even better. I found three other tools:

1. MuSiC: https://www.nature.com/articles/s41467-018-08023-x
Published in nature communication and a few months older than cibersortx. I haven't studied their algorithm yet and so I can't say if their method is better than cibersortx or not. However, the part that makes the difference is that all of their works, including their source code, is open-source. So, I can access their code and we don't need to upload our unpublished data into some not safe online server (unlike cibersortx). The second good thing is that it is implemented in R language which I use for most of my RNA-seq data analysis, so it is easily compatible with other tools.
2. Scaden : https://advances.sciencemag.org/content/6/30/eaba2619.full This one is also open-source and accessible easily. It is written in Python and runs in the shell. So, it is somewhat compatible with other tools (not as compatible as MuSic). We again can process our data locally and don't need to upload our data to any servers (which is good). They claim that their method has a better correlation with the real scRNA-seq data than both cibersortx and MuSic (not always!!). Running their tools seems easy based on their tutorials (I haven't run it on my system yet!). However, their tutorials seem to be a little confusing and it seems like it is still in pre-release version and they have bugs!!
3. CDSeq: https://github.com/kkang7/CDSeq\_R\_Package Finally, I found the CDSeq package which doesn't need any scRNA-seq input to find single cell proportions and cell-specific gene expression in bulk data. It just gets your buld counts as input and info about very few parameters and does the whole thing automatically, and suprisingly good, I checked! If you could find gene expression pattern of your specific cells and also use it as an input, it is optional, it will give you even a better estimation.
> 0 votes

---
Tags: rna-seq, cell-line
---


# Question
Title: NO_COOR reads not in a single block at the end 0 -1
I am doing an RNA seq analysis with mouse samples. I had my reads in several lanes, so I concatenated those fastq reads. Then I've mapped my reads with Hisat2 against the new version of mouse "GCA\_000001635.9\_GRCm39\_genomic.fna".

Apparently everything works well up to here. Now I'm using samtools 1.9 to convert sam to bam, to sort those bam and create the index, just like this:

```
samtools view -bS ${file} > ${file}.bam

samtools sort -n ${file} -o ${file}_sorted.bam

samtools index  ${file}_sorted.bam

```

But when indexing I get this error:

> \[bam\_sort\_core\] merging from 11 files and 1 in-memory blocks... \[E::hts\_idx\_push\] NO\_COOR reads not in a single block at the end 0 -1 samtools index: failed to create index for "CD4\_1A\_concat.bam\_sorted.bam": No such file or directory
>
> \[bam\_sort\_core\] merging from 11 files and 1 in-memory blocks... \[E::hts\_idx\_push\] NO\_COOR reads not in a single block at the end 23 -1 samtools index: failed to create index for "CD4\_1B\_concat.bam\_sorted.bam": No such file or directory \[bam\_sort\_core\] merging from 11 files and 1 in-memory blocks...

And more...

Any help would be appreciated!

# Answer
This is essentially a duplicate of Samtools Index: Chromosome Blocks not Continuous but with different error messages, and a duplicate of Error given while trying to index a BAM file with Samtools Index - NO COOR?.

```
samtools sort -n ${file} -o ${file}_sorted.bam

```

These error messages also indicate that the reads are not sorted by coordinate — in particular, that the unplaced unmapped reads with RNAME `*` are not all together at the end of the file.

This is because `samtools sort -n` has been used to sort the reads by name instead. Remove `-n` to sort by position, which is what is needed to prepare a BAM file for indexing with `samtools index`.
> 2 votes

---
Tags: rna-seq, samtools
---


# Question
Title: Master of Pores Installation Mac OS Mojave 10.14
I'm bit of a newbie with ONT analyses and am having some trouble with the Master Of Pores installation instructions (https://biocorecrg.github.io/master\_of\_pores/install.html). I've managed to install nextflow, cloned the repository and have downloaded albacore and guppy. But when I run the pipeline below, i get an error saying nextflow: command not found. Any help would be most appreciated.

```
curl https://get.nextflow.io | bash
mv nextflow ~/master_of_pores/NanoPreProcess/bin
nextflow run nanopreprocess.nf

```

# Answer
> 3 votes
you should not move nextflow within `~/master_of_pores/NanoPreProcess/bin` . Just add nextflow in your path.

```
export PATH=$PATH:where_is_your_nf_executable

```

Feel free to open an issue on github if you need more help.

---
Tags: nanopore, software-installation, nextflow
---


# Question
Title: Tracing and quantifying inheritance over many generations
Apologies if this is the wrong location for this question, and please feel free to point me in the right direction with appropriate Anglo-Saxon embellishments...

I'm after some ball-park numbers for the following questions (plus any hints how to calculate them for different params - e.g. less generations, larger starting pop.):

Given a starting population of 10,000, and given 2,500 generations, with an end population of 5,000,000,000

1. What are the chances of being a direct descendant of 1 individual of the original 10,000?
2. If you had a DNA sample of that individual, would you be able to tell if you were a direct descendant?
3. Would siblings of the original individual confuse things with regard to Q2?

# Answer
> 0 votes
Since it seems my comment helped, I've been suggested to post this as an answer...

I'd reccomend everyone interested in these sort of population genetics to have a look at 'The Ancestors Tale' by Richard Dawkins. In particular, chapter 0, which deals with human/hominid genetics. That should help answer questions. The relevant part of the book appears to be available online for free in GoogleBooks (https://books.google.co.uk/books/about/The\_Ancestor\_s\_Tale.html?id=s86\_CQAAQBAJ).

Below I'll note down the answers from the comments section, not because I wrote them but I figured it might help anyone who stumbled over this question afterwards:

1. The chances of being a direct descendant of 1 individual of the original 10,000 is 100%, assuming they have any descendants. - Tom Melly
2. If one had a DNA sample of that individual, one could use uniparental markers (either Y chromosome for men, or mtDNA for women) to infer whether they are a direct descendant. - 4galaxy7

(My guess is that it doesn't have to be uniparental markers. Other genes might work just as well. I'm thinking of Neanderthal/Denisovan work. But I might be wrong.)

# Answer
> 0 votes
## Question 1

If I understand the question, then you are asking 'Given a starting population of 10 000 individuals and 2 500 generations, what is the probability that one of the 5 000 000 000 individuals descendend individuals is directly descended from one of the founders?'.

If we are assuming no migration in or out of the system, then the answer must be 1, since there are no ways to be descended from someone outside the system. Please let me know if I have misunderstood the question.

## Question 2

For this, it's only possible to use uni-parental markers. This is because inferring exact genealogical relationships using autosomal markers beyond something like 3rd cousins is extremely difficult, owing to the stochastic way in which recombination breaks up segments. Beyond a few generations, only the uniparentally inherited mtDNA and non-recombining Y-chromosome can be informative about relatedness.

Uniparental markers don't recombine, so mutations aside, you expect to inherit an exact copy of your mothers mtDNA marker. If you know the number of generations which are likely to separate you and the putative ancestor, you can estimate the some kind of likelihood of being directly descended from the ancestor given the number of pairwise differences between the individual and the ancestor.

This is what they did when they wanted to identify the remains of Richard III. It's a good paper and easy to read / open access, so I'd recommend checking it out.

## Question 3

Simply, no. The inference of being a direct descendent from one of the ancestors would be independent from any of the other potential descendants in the sample.

---
Tags: phylogenetics
---


# Question
Title: TCGA: gene IDs not appearing and other concerns
I want to download RNA-seq datasets of 4 or 5 different types of cancer from the TCGA to investigate what is happening with my gene of interest. The problem is that I'm a physician and I don't know people used to deal with Big Data to help me.

I'm using a Bioconductor program called TCGAbiolinks and I'm able to download some datasets using the following script:

```
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install("TCGAbiolinks")

browseVignettes("TCGAbiolinks")

library(TCGAbiolinks)
library(dplyr)
library(DT)

query <- GDCquery(project = "TCGA-BRCA",
                  data.category = "Gene expression",
                  file.type  = "normalized_results",
                  data.type = "Gene expression quantification",
                  experimental.strategy = "RNA-Seq",
                  barcode = c("TCGA-AN-A046", "TCGA-AC-A23H"),
                  legacy = TRUE)

GDCdownload(query, method = "api", files.per.chunk = 10)

# GDCprepare transforms data in R format

# The default GDCprepare output is a SummarizedExperiment (SE) object,
# but you can also set to output the data frame.

data <- GDCprepare(query, summarizedExperiment = FALSE)

```

The problems is that: 1) the gene list is not appearing in the first column; 2) I don't know how many datasets I'll be able to download and analyse without a server.. I'm working with my Macbook..

Can someone familiar with the TCGA website or with the TCGAbiolinks help me with this project? I've been trying to progress in the last months, but the TCGA is a complex website and I'm feeling lost. I use R, I don't know SQL or other programming languages..

# Answer
I had to play around a little to understand the object, but here's what I got:

```
R version 4.0.2 (2020-06-22) -- "Taking Off Again"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin17.0 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(TCGAbiolinks)
> query <- GDCquery(project = "TCGA-BRCA",
+                   data.category = "Gene expression",
+                   file.type  = "normalized_results",
+                   data.type = "Gene expression quantification",
+                   experimental.strategy = "RNA-Seq",
+                   barcode = c("TCGA-AN-A046", "TCGA-AC-A23H"),
+                   legacy = TRUE)
--------------------------------------
o GDCquery: Searching in GDC database
--------------------------------------
Genome of reference: hg19
--------------------------------------------
oo Accessing GDC. This might take a while...
--------------------------------------------
ooo Project: TCGA-BRCA
--------------------
oo Filtering results
--------------------
ooo By experimental.strategy
ooo By data.type
ooo By file.type
ooo By barcode
----------------
oo Checking data
----------------
ooo Check if there are duplicated cases
ooo Check if there results for the query
-------------------
o Preparing output
-------------------

```

---

```
> GDCdownload(query, method = "api", files.per.chunk = 10)
Downloading data for project TCGA-BRCA
GDCdownload will download 3 files. A total of 1.306357 MB
Downloading chunk 1 of 1 (3 files, size = 1.306357 MB) as Wed_Aug_19_12_20_11_2020_0.tar.gz

Downloading: 580 kB

```

---

```
> data <- GDCprepare(query, summarizedExperiment = TRUE)
------------------
oo Reading 3 files
------------------
|====================================================|100%                      Completed after 0 s
------------------
oo Merging 3 files
------------------
Accessing grch37.ensembl.org to get gene information
Downloading genome information (try:0) Using: Human genes (GRCh37.p13)
Starting to add information to samples
 => Add clinical information to samples
 => Adding TCGA molecular information from marker papers
 => Information will have prefix 'paper_'
brca subtype information from:doi.org/10.1016/j.ccell.2018.03.014

```

---

```
> head(data)
class: RangedSummarizedExperiment
dim: 6 3
metadata(1): data_release
assays(1): normalized_count
rownames(6): A1BG A2M ... RP11-986E7.7 AADAC
rowData names(3): gene_id entrezgene ensembl_gene_id
colnames(3): TCGA-AN-A046-01A-21R-A034-07 TCGA-AC-A23H-11A-12R-A157-07
  TCGA-AC-A23H-01A-11R-A157-07
colData names(83): barcode patient ... paper_PARADIGM Clusters
  paper_Pan-Gyn Clusters

```

I see from this representation that I can use `metadata(data)`, `assays(data)`, etc to get more information.

```
> head(SummarizedExperiment::assays(data))
List of length 1
names(1): normalized_count

```

Alright, so this is a named list. I know how to deal with lists: I can use either the index (`[[1]]`) or the name (`$normalized_count`)

```
> head(SummarizedExperiment::assays(data)$normalized_count)
             TCGA-AN-A046-01A-21R-A034-07 TCGA-AC-A23H-11A-12R-A157-07 TCGA-AC-A23H-01A-11R-A157-07
A1BG                             121.9879                      38.5795                      13.6874
A2M                             9703.8622                   37640.8373                    4767.0490
NAT1                            8488.0727                     114.8582                     270.4753
NAT2                               6.0583                       4.2938                       2.3804
RP11-986E7.7                   33016.6641                   11308.3460                    3894.9639
AADAC                              1.1359                     219.6976                       1.1902

```

I can now get more information on the samples using `colData()`. I'm truncating the output below because there are 83 columns.

```
> SummarizedExperiment::colData(data)
DataFrame with 3 rows and 83 columns
                                                  barcode      patient           sample shortLetterCode          definition sample_submitter_id         sample_type       state
                                              <character>  <character>      <character>     <character>         <character>         <character>         <character> <character>
TCGA-AN-A046-01A-21R-A034-07 TCGA-AN-A046-01A-21R-A034-07 TCGA-AN-A046 TCGA-AN-A046-01A              TP Primary solid Tumor    TCGA-AN-A046-01A       Primary Tumor    released
TCGA-AC-A23H-11A-12R-A157-07 TCGA-AC-A23H-11A-12R-A157-07 TCGA-AC-A23H TCGA-AC-A23H-11A              NT Solid Tissue Normal    TCGA-AC-A23H-11A Solid Tissue Normal    released
TCGA-AC-A23H-01A-11R-A157-07 TCGA-AC-A23H-01A-11R-A157-07 TCGA-AC-A23H TCGA-AC-A23H-01A              TP Primary solid Tumor    TCGA-AC-A23H-01A       Primary Tumor    released
                                                        sample_id   is_ffpe sample_type_id  tissue_type submitter_id days_to_collection oct_embedded                pathology_report_uuid
                                                      <character> <logical>    <character>  <character>  <character>          <integer>  <character>                          <character>
TCGA-AN-A046-01A-21R-A034-07 57aff09f-0e97-4e65-a227-3dc7a8516367     FALSE             01 Not Reported TCGA-AN-A046                 27         true 1304FB17-A20A-4EBC-9CC8-1554808AC1F6
TCGA-AC-A23H-11A-12R-A157-07 7df59ca8-7e51-4581-9d7f-8bba0395ce17     FALSE             11 Not Reported TCGA-AC-A23H                478        false                                   NA
TCGA-AC-A23H-01A-11R-A157-07 d7e3b628-d5fd-4e79-9c4a-6409330fb8a7     FALSE             01 Not Reported TCGA-AC-A23H                478        false A7C7D409-D086-4A9B-8C8F-E7E231D5891D

```

I can get more information on the genes using `rowData()`:

```
> SummarizedExperiment::rowData(data)
DataFrame with 19947 rows and 3 columns
                  gene_id entrezgene ensembl_gene_id
              <character>  <integer>     <character>
A1BG                 A1BG          1 ENSG00000121410
A2M                   A2M          2 ENSG00000175899
NAT1                 NAT1          9 ENSG00000171428
NAT2                 NAT2         10 ENSG00000156006
RP11-986E7.7 RP11-986E7.7         12 ENSG00000273259
...                   ...        ...             ...
RASAL2-AS1     RASAL2-AS1  100302401 ENSG00000224687
LINC00882       LINC00882  100302640 ENSG00000242759
FTX                   FTX  100302692 ENSG00000230590
TICAM2             TICAM2  100302736 ENSG00000243414
SLC25A5-AS1   SLC25A5-AS1  100303728 ENSG00000224281

```

You want the assay data as a data.matrix, I guess. That'd be:

```
> data_as_matrix <- data.matrix(SummarizedExperiment::assays(data)$normalized_count)

```

My session Info:

```
> sessionInfo()
R version 4.0.2 (2020-06-22)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.5

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] TCGAbiolinks_2.16.1

loaded via a namespace (and not attached):
  [1] bitops_1.0-6                matrixStats_0.56.0          bit64_0.9-7.1               doParallel_1.0.15           RColorBrewer_1.1-2          progress_1.2.2
  [7] httr_1.4.2                  GenomeInfoDb_1.24.2         backports_1.1.8             tools_4.0.2                 R6_2.4.1                    DBI_1.1.0
 [13] BiocGenerics_0.34.0         colorspace_1.4-1            tidyselect_1.1.0            gridExtra_2.3               prettyunits_1.1.1           bit_1.1-15.2
 [19] curl_4.3                    compiler_4.0.2              rvest_0.3.6                 Biobase_2.48.0              xml2_1.3.2                  DelayedArray_0.14.0
 [25] rtracklayer_1.48.0          scales_1.1.1                survMisc_0.5.5              readr_1.3.1                 genefilter_1.70.0           askpass_1.1
 [31] rappdirs_0.3.1              stringr_1.4.0               digest_0.6.25               Rsamtools_2.4.0             foreign_0.8-80              R.utils_2.9.2
 [37] rio_0.5.16                  XVector_0.28.0              pkgconfig_2.0.3             dbplyr_1.4.4                readxl_1.3.1                rlang_0.4.7
 [43] ggthemes_4.2.0              RSQLite_2.2.0               generics_0.0.2              zoo_1.8-8                   jsonlite_1.7.0              BiocParallel_1.22.0
 [49] zip_2.0.4                   car_3.0-8                   dplyr_1.0.0                 R.oo_1.23.0                 RCurl_1.98-1.2              magrittr_1.5
 [55] GenomeInfoDbData_1.2.3      Matrix_1.2-18               Rcpp_1.0.5                  munsell_0.5.0               S4Vectors_0.26.1            abind_1.4-5
 [61] lifecycle_0.2.0             R.methodsS3_1.8.0           stringi_1.4.6               carData_3.0-4               SummarizedExperiment_1.18.1 zlibbioc_1.34.0
 [67] plyr_1.8.6                  BiocFileCache_1.12.0        grid_4.0.2                  blob_1.2.1                  parallel_4.0.2              ggrepel_0.8.2
 [73] forcats_0.5.0               crayon_1.3.4                survminer_0.4.8             lattice_0.20-41             haven_2.3.1                 Biostrings_2.56.0
 [79] splines_4.0.2               GenomicFeatures_1.40.0      annotate_1.66.0             hms_0.5.3                   knitr_1.29                  pillar_1.4.6
 [85] ggpubr_0.4.0                GenomicRanges_1.40.0        ggsignif_0.6.0              codetools_0.2-16            biomaRt_2.44.1              stats4_4.0.2
 [91] XML_3.99-0.5                glue_1.4.1                  downloader_0.4              data.table_1.13.0           vctrs_0.3.2                 selectr_0.4-2
 [97] foreach_1.5.0               cellranger_1.1.0            gtable_0.3.0                openssl_1.4.2               purrr_0.3.4                 tidyr_1.1.0
[103] km.ci_0.5-2                 assertthat_0.2.1            ggplot2_3.3.2               openxlsx_4.1.5              xfun_0.16                   xtable_1.8-4
[109] broom_0.7.0                 rstatix_0.6.0               survival_3.1-12             tibble_3.0.3                iterators_1.0.12            KMsurv_0.1-5
[115] GenomicAlignments_1.24.0    AnnotationDbi_1.50.1        memoise_1.1.0               IRanges_2.22.2              ellipsis_0.3.1

```
> 3 votes

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