title: BLASTn query coverage discrepancy

text:
# Question
Title: BLASTn query coverage discrepancy
Using the BLASTn webtool with,

**Query Seq:** CAATTGCGAAGCCACATTATT

**Subject Seq:** CP031332.1

BLASTn returns a single hit to CP031332.1 with query coverage at 100%.

However if one looks at the alignment, the query coverage is clearly 17/21 ≈ 81%. Could someone explain to me the reason for this discrepancy?

Edit,

Web version: https://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastn&PAGE\_TYPE=BlastSearch&LINK\_LOC=blasthome

Output: https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Get&RID=RCK3F55Y114

# Answer
> 1 votes
To expand on the excellent answer by *haci*:
BLAST reports (a) the summary of alignments and (b) individual alignments, see below. Note that your link to the results expired (as NCBI BLAST results links tend to do). It is considered best practice to post the text version of the relevant portion of the results in the question itself.

The summary provides an overview of all the query sequence **individual hits** in the subject sequence. Because the individual hits (alignments, or ranges) **together** cover 100% of the query and each of these has the max identity = 100%, that is what it reports in the summary. It is not very intuitive for many users.

You should always look at the individual hits in detail. The above 100% coverage in the summary consists of ranges 1 and 14, shown below. The top hit (range 1) covers only 17 out of 21 nt with 100% identity, and this is the most relevant information from BLAST. The other hit (range 14) covers 10 out of 21 nt with 100% identity. This is the cause of the discrepancy.

These 2 alignments are not close together in the subject sequence. Depending on your goal, choose one or more of the individual alignments and consider alignment length and % identity of those hits individually, rather than summarized together as in the BLAST report summary. For most applications, the top hit (range 1) would be the most relevant alignment.

```
Program: BLASTN
Subject:Neisseria meningitidis strain M22814 chromosome, complete genome ID: CP031332.1(nucleic acid) Length: 2190201
Query #1: seq1 Query ID: lcl|Query_56011 Length: 21

Sequences producing significant alignments:
                                                                  Max    Total Query   E   Per.
Description                                                       Score  Score cover Value Ident  Accession
Neisseria meningitidis strain M22814 chromosome, complete genome  34.2   57729 100%  0.001 100.00 CP031332.1

Alignments:
>Neisseria meningitidis strain M22814 chromosome, complete genome
Sequence ID: CP031332.1 Length: 2190201
Range 1: 2074883 to 2074899

Score:34.2 bits(17), Expect:0.001,
Identities:17/17(100%),  Gaps:0/17(0%), Strand: Plus/Plus

Query  1        CAATTGCGAAGCCACAT  17
                |||||||||||||||||
Sbjct  2074883  CAATTGCGAAGCCACAT  2074899

...

Range 14: 354292 to 354301

Score:20.3 bits(10), Expect:15,
Identities:10/10(100%),  Gaps:0/10(0%), Strand: Plus/Minus

Query  12      CCACATTATT  21
               ||||||||||
Sbjct  354301  CCACATTATT  354292

```

# Answer
> 0 votes
Your blast result is not composed of a single alignment of all of the length of your query vs the subject but from a series of alignments: in each a subsequence of your query is used and these alignments are sorted by "E value", "score", "percent identity", ...

Now in your case, you are talking about the first hit, in which a subsequence of your query (positions 1-17) match perfectly to a portion of the subject. Query coverage is not calculated from this first hit though, it takes all of the alignments into account and again in your case all 21 bases were used in one or more alignments. I believe it would be fair to filter out all but the first hit given the E-values and re-calculate the query coverage as 17/21.

---
Tags: blast, blastn
---


title: Easy way to create gene plot in gviz?

text:
# Question
Title: Easy way to create gene plot in gviz?
I would like to create a simple plot of a given human gene containing the main features, e.g. exons, introns, ORFs, transcripts etc. I understand that `gviz` is a popular choice for this kind of thing, but most tutorials are massive overkill for what I want to do (I don't need fancy features such as importing BAM files and plotting coverage).

I was wondering whether someone can point me towards some sample code for creating a basic gene plot. It probably would be good to start by creating a dataframe containing all the features of my gene - is there a handy package that can do that?

# Answer
> 1 votes
If Python is an option, you could check BioPython and GenomDiagram. This should get the job done.

# Answer
> 0 votes
I assume you are trying to make something like a "gene arrow map" and `[gggenes][1]` can help you (with some little tweaks) in case you are familiar with `ggplot2`.

---
Tags: gene, visualization, gviz
---


# Question
Title: Error while using E-utilities on slurm HPC
I keep getting the error while using E-utilities on Slurm HPC, could you suggest how I can resolve it?

SCRIPT:

```
#!/bin/bash
#SBATCH --job-name=tax_job_test      # Job name
#SBATCH --mail-type=END,FAIL         # Mail events (NONE, BEGIN, END, FAIL, ALL)
#SBATCH --mail-user=na@gmail.com    # Where to send mail
#SBATCH --ntasks=8                  # Number of MPI ranks
#SBATCH --nodes=1                    # Number of nodes
#SBATCH --time=23:05:00              # Time limit hrs:min:sec
#SBATCH --partition=k2-medpri

module load e-utilities/03.02.19
module load apps/perl/5.20.2/gcc-4.8.5

for i in `cat taxonomy.id.txt` ; do efetch -db nuccore -format gbc -id ${i} |xtract  -pattern INSDSeq -element INSDSeq_accession-version INSDSeq_taxonomy INSDSeq_organism >>taxo.2020.txt; done

```

input `taxonomy.id.txt` is like this:

```
AKN43077.1
AKN41814.1
PNE10086.1
AKN42672.1
AKN41654.2
AKN43077.1
PNE10086.1
AKN43189.1
TMK08501.1
TMK08501.1
AKN42976.1
AKN43189.1
AIZ47938.1

```

Error:

The script keep giving below error:

```
501 Protocol scheme 'https' is not supported (LWP::Protocol::https not installed)
No do_post output returned from 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=AKN43077.1
&rettype=gbc&retmode=xml&edirect_os=linux&edirect=11.7&tool=edirect&email=3052771@node126.pri.kelvin2.alces.network'

Result of do_post http request is

$VAR1 = bless( {
                 '_rc' => 501,
                 '_request' => bless( {
                                        '_method' => 'POST',
                                        '_uri' => bless( do{\(my $o = 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi')}, 'URI::https' ),
                                        '_content' => 'db=nuccore&id=AKN43077.1
&rettype=gbc&retmode=xml&edirect_os=linux&edirect=11.7&tool=edirect&email=3052771@node126.pri.kelvin2.alces.network',
                                        '_headers' => bless( {
                                                               'content-type' => 'application/x-www-form-urlencoded',
                                                               'user-agent' => 'libwww-perl/6.39'
                                                             }, 'HTTP::Headers' )
                                      }, 'HTTP::Request' ),
                 '_content' => 'LWP will support https URLs if the LWP::Protocol::https module
is installed.
',
                 '_headers' => bless( {
                                        'client-warning' => 'Internal response',
                                        '::std_case' => {
                                                          'client-warning' => 'Client-Warning',
                                                          'client-date' => 'Client-Date'
                                                        },
                                        'client-date' => 'Mon, 05 Oct 2020 14:27:37 GMT',
                                        'content-type' => 'text/plain'
                                      }, 'HTTP::Headers' ),

```

Many thanks

# Answer
> 2 votes
`501 Protocol scheme 'https' is not supported (LWP::Protocol::https not installed)` indicates that a specific perl module is missing.

If indeed your edirect installation is from March 2019 (as I interpret from `e-utilities/03.02.19`), you should update your edirect installation to the latest version. If you update, I expect the necessary perl module will be installed and the error will go away.

---
Tags: phylogenetics, hpc, e-utilities
---


# Question
Title: Bootstrapping cell type identifications (scRNA-seq)
I have scRNA-seq data held in a `SEURAT` object, and cell-types were identified with `SingleR` (essentially a vector of strings - Monocytes, Granulocytes etc.). I have two conditions - treated and untreated. I want to know if the difference in the relative frequency (%) of a given cell type between the conditions is significant. For this, I thought to perform bootstrapping, to obtain frequencies for many different 'simulated' datasets - and then perform a statistical test. However, I cannot understand how best to do this. Should I subsample the scRNA-seq data, with replacement and calculate the frequencies each time? Can somebody provide a step by step workflow on how to do this?

# Answer
I would use permutation:

i) Get the difference that you are interested in, i.e. number/proportion of monocytes in one group vs the other.

ii) Get the relevant columns (condition, cell label, ...) from the `@meta.data` slot of your seurat objects and `rbind()`.

iii) Shuffle the cell labels randomly and split the data frame into two (per condition).

iv) Compare the difference between the cell counts/proportions of interest from the two permuted data frames and save this value.

iv) Repeat iii and iv for 100s of times and calculate how many times have you observed a difference at least as extreme as your real observation (i). Once divided by the number of trials, this would give you a p-value stating how likely it is that you get the observed value.

I have to add that batch effects are severe in scRNA-seq and hence a linear model where you might be able to account for batches could be more helpful.
> 0 votes

---
Tags: scrnaseq, statistics
---


# Question
Title: How to identify unknow bacterium species from whole genome genetic sequence
I am a Biochemist that is unfamiliar with bioinformatic tools and new to academia as a whole.

I am currently using ILLUMINA PE data, which I trimmed (Trimmomatic), corrected (Rcorrector) and assembled (SPAdes). I am now interested in using the genetic sequences from my contigs to identify the original organism (I know it is a bacterium).

I have tried using BLASTn (with default parameters) on one of my contigs. I let the query run for 30 minutes before stopping, nothing came of it.

What would be the best way to go about doing this?

# Answer
You can use KmerFinder and Kraken2 for taxonomy assignment
> 1 votes

---
Tags: fasta, blast, genomics, blastn, genetics
---


# Question
Title: Seurat Violin Plot: Why do dots align in one row?
I am analyzing a rather small scRNAseq dataset (39797 features across 3522 samples) and I have observed that dots in the violin plot are plotted in one line next to each other rather than in a cloud of dots (see plot below):

I am just wonndering what the reason might be for this?

I have followed the Integration and Label Transfer tutorial using SCTransform (https://satijalab.org/seurat/v3.1/integration.html). After the integration process I switched back to the SCT assay, which should contain the normalized counts. For plotting I used the following command: VlnPlot(..., assay = "SCT", slot = "data"), again specifying to use the normalized counts.

Can someone explain me why this is the case?

Thank you already in advance!

# Answer
That is because SCTransform (mostly) outputs discrete values, you can check the @slot of your Seurat object.
> 1 votes

---
Tags: scrnaseq, seurat
---


title: How many entries in PDB have SS bonds?

text:
# Question
Title: How many entries in PDB have SS bonds?
Is there a way to figure out what percentage of entries in PDB have SS bonds? In addition, what other types of bonds in PDB may be worth looking into (not H bonds)?

# Answer
The easiest way is to check annotations in mmCIF or PDB files.

Looking at my page with statistics of mmCIF tags, last updated in March:

```
_struct_conn_type.id        59137×metalc 42706×covale 36195×disulf 10725×hydrog
_struct_conn.conn_type_id   5350402×hydrog 1733277×metalc 559092×covale 226472×disulf

```

Which means 226472 disulfide bonds in 36195 PDB entries at the time when the PDB had 161470 entries.
> 3 votes

---
Tags: proteins, database, protein-structure, pdb
---


title: alpha diversity wilcox.test

text:
# Question
Title: alpha diversity wilcox.test
I am trying to do Wilcox test to detect the significant difference in alpha diversity but it is showing an error? physeqN2 is a phyloseq object and Season is a metadata column.

```
results = estimate_richness(physeqN2, measures = 'Chao1')

d = sample_data(physeqN2)

# calculate Wilcox Test
Wet.root = results[d[,'Season'] == 'Wet',]
Dry.root = results[d[,'Season'] == 'Dry',]
pmv <-pairwise.wilcox.test(Wet.root, Dry.root)

capture.output(pmv, file = "chao1.16s.wet_vs_Dry.root.txt")

Error;
Error in wilcox.test.default(xi, xj, paired = paired, ...) :
  'x' must be numeric
Calls: pairwise.wilcox.test ... compare.levels -> wilcox.test -> wilcox.test.def                                                                                ault
Execution halted

```

I hvae tried to do it in another way but here also getting error?

```
alpha.diversity <- estimate_richness(physeqN2, measures=c("Chao1"))

## Wilcox on observed richness
data <- cbind(sample_data(physeqN2), alpha.diversity)
pmv <-pairwise.wilcox.test(Chao1 ~ Season, data)

capture.output(pmv, file = "chao1.16s.wet_vs_Dry.root.txt")

Error:

Error in if (length(ans) == 0L || as.character(ans[[1L]])[1L] == "~") { :
  missing value where TRUE/FALSE needed
Calls: pairwise.wilcox.test ... sapply -> lapply -> FUN -> compare.levels -> [ -> [.formula
Execution halted

```

Many thanks

# Answer
I think this is a bug in the input data.

At a wild guess what I think is happening is the input data is not regarded as paired, hence the error `paired = paired`. The other weird error is `x must be numeric`. Firstly,

* how are you scoring missing data, its not via `-` by any chance? Missing data must be via `nan` type syntax.
* are you sure there are no heading bugs? Like a header is getting dragged into the numeric data?

Wilcoxan is a non-parametric T test, so the wet root and dry root must match exactly the same sample. Thus it is not a general wet root versus generic dry root, but wet and dry from the same specimen. If this is the data you have the program doesn't think it is matching.

If it is truely paired, again I think the source of error is you have some missing value which is leaving the pairings unequal lengths or else a header bug.

If the data genuinely isn't paired, then I would look at the Mann-Whitney U test, it would be easier in this scenario.

I don't use R syntax (just pandas dataframes) so can't say whether the syntax in the first example is valid, but it is pretty clear regarding your input. In Python you'd generally just give it the entire dataframe and tell it to go figure.

The reality is that loads of programs do Wilcoxan, you don't need to use R if its a one-off test. Obviously if this a part of a pipeline (automation) thats a bit different.
> 0 votes

---
Tags: r, bioconductor, statistics, microbiology
---


# Question
Title: data visualization RNAseq : scaling data for PCA and cluster dendogram
I have count data from a RNAseq experiment (2 samples are from normal cells and 3 samples are cells with a disease), and the data is already standardized by trimmed mean of M values (TMM). I want to do some plots: biplot of Principal Component Analysis (PCA) and a cluster dendrogram to see if the samples normal/disease are well separated (there is a clear difference between them). Since the data is already standardized (TMM) should I scale and center the data prior to perform PCA and cluster dendrogram??

thank you!

# Answer
> 1 votes
Scaling (or centering) makes the genes comparable: Putting the expression levels of genes in the same scale (i.e. between 0 and 1) sustains that all of your genes contribute equally to the PCA or distance calculations. On the other hand, without this step, such calculations would be dominated by the highly/lowly expressed genes. I get "better" results with scaling.

---
Tags: rna-seq, clustering, pca, tmm
---


# Question
Title: gene-gene correlation from two different Tissues
I am stuck on how to do correlation for two independent data sets with common row and column names. A and B are datasets that contain as many rows as genes and as many columns as samples. The rows in A and B represent a common set of genes but measured in two different tissues. The columns represent measurements in the same 5 samples in both A and B. I want to do a correlation between the set of genes in A and B. This is to see if the same genes in both tissues are correlated or not. Since the matrix would be big in my actual data, I only want to retain a correlation coefficient higher than 0.5.

Here I simulate the data set.

```
set.seed(1)
A <- data.frame(rnorm(100),
                rnorm(100),
                rnorm(100),
                rnorm(100),
                rnorm(100))
row.names(A) <- paste0("G_", 1:100)
colnames(A) <- paste0("M_", 1:5)

set.seed(42)
B <- data.frame(rnorm(100),
                rnorm(100),
                rnorm(100),
                rnorm(100),
                rnorm(100))
row.names(B) <- paste0("G_", 1:100)
colnames(B) <- paste0("I_", 1:5)

```

Thank you!

# Answer
> 1 votes
You can use `mapply()`. `ta` and `tb` being transposed data frames of your `A` and `B` data frames respectively:

```
> mapply(cor, ta, tb)[mapply(cor, ta, tb) > 0.5]

      G_3       G_5       G_9      G_10      G_11      G_15      G_20      G_23      G_25      G_26      G_33      G_40      G_43      G_48      G_57      G_60
0.5346591 0.8066507 0.8379777 0.6752681 0.7221359 0.5285787 0.7333045 0.5627962 0.6533379 0.7256878 0.5996492 0.6486557 0.5108215 0.7386332 0.6596823 0.6919915
     G_63      G_72      G_76      G_80      G_81      G_90      G_97      G_98      G_99
0.5589583 0.8391917 0.7608801 0.8003665 0.6364557 0.5030968 0.7298439 0.5693024 0.5709411

```

---
Tags: r, gene, gene-expression, correlation
---


# Question
Title: How do I obtain a vcf read from a pyvcf.Reader?
I wondered how would I be able to obtain variant based on chrom and location range.

Optimally, this function would be like:

```
import vcf
vcf_reader = vcf.Reader(filename=some_file_name)
vcf_read = vcf.variants_between(CHROM, START, END)

```

Does anyone how to do so?

# Answer
> 1 votes
I think `vcf_reader.fetch` is what you want. You have to index the file first though.

If starting from regular text VCF you can run:

```
bgzip example.vcf
tabix example.vcf.gz

```

Then the Python code is almost the same as yours:

```
import vcf
vcf_reader = vcf.Reader(filename="example.vcf.gz")
vcf_read = vcf_reader.fetch("NC_023663.1", 7803980, 7804000)
for record in vcf_read:
    print(record)

```

For this example VCF you get:

```
Record(CHROM=NC_023663.1, POS=7803981, REF=G, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803982, REF=T, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803983, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803984, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803985, REF=G, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803986, REF=G, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803987, REF=G, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803988, REF=C, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803989, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803990, REF=T, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803991, REF=A, ALT=[G, <*>])
Record(CHROM=NC_023663.1, POS=7803992, REF=G, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803993, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803994, REF=C, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803995, REF=T, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803996, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803997, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803998, REF=C, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803999, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7804000, REF=T, ALT=[<*>])

```

API reference: https://pyvcf.readthedocs.io/en/latest/API.html#vcf.Reader.fetch

Note the "zero-based, half-open coordinate system" it describes (which is why 7803980 in the second Python argument becomes POS=7803981 from the VCF but then it stops at POS=7804000).

---
Tags: vcf, python
---


title: amino acids substitution probability

text:
# Question
Title: amino acids substitution probability
I am performing amino acid analysis but require assistance in the underlying amino acid substitution model, i.e. the predetermined rate/frequency of amino acid mutations which are used for example in constructing phylogenetics trees using maximum likelihood or Bayesian methods.

There are numerous matrices of probabilities for every combination of amino acids (something like *BLOSUM* or *PAM* matrices and I am wondering within the large number of models which the most appropriate. I am aware there is no 'universal' matrix.

I want to use the matrix to check the probability of amino acid substitution in sequence on a pairwise basis.

For example, I'm not sure whether the BLOSUM or PAM matrix will be useful to me.

\[**EXTRA**\] If all diagonal values in such a matrix would be equal to `1` that would be helpful.

# Answer
The LG model is here:

```
0.425093
0.276818 0.751878
0.395144 0.123954 5.076149
2.489084 0.534551 0.528768 0.062556
0.969894 2.807908 1.695752 0.523386 0.084808
1.038545 0.363970 0.541712 5.243870 0.003499 4.128591
2.066040 0.390192 1.437645 0.844926 0.569265 0.267959 0.348847
0.358858 2.426601 4.509238 0.927114 0.640543 4.813505 0.423881 0.311484
0.149830 0.126991 0.191503 0.010690 0.320627 0.072854 0.044265 0.008705 0.108882
0.395337 0.301848 0.068427 0.015076 0.594007 0.582457 0.069673 0.044261 0.366317 4.145067
0.536518 6.326067 2.145078 0.282959 0.013266 3.234294 1.807177 0.296636 0.697264 0.159069 0.137500
1.124035 0.484133 0.371004 0.025548 0.893680 1.672569 0.173735 0.139538 0.442472 4.273607 6.312358 0.656604
0.253701 0.052722 0.089525 0.017416 1.105251 0.035855 0.018811 0.089586 0.682139 1.112727 2.592692 0.023918 1.798853
1.177651 0.332533 0.161787 0.394456 0.075382 0.624294 0.419409 0.196961 0.508851 0.078281 0.249060 0.390322 0.099849 0.094464
4.727182 0.858151 4.008358 1.240275 2.784478 1.223828 0.611973 1.739990 0.990012 0.064105 0.182287 0.748683 0.346960 0.361819 1.338132
2.139501 0.578987 2.000679 0.425860 1.143480 1.080136 0.604545 0.129836 0.584262 1.033739 0.302936 1.136863 2.020366 0.165001 0.571468 6.472279
0.180717 0.593607 0.045376 0.029890 0.670128 0.236199 0.077852 0.268491 0.597054 0.111660 0.619632 0.049906 0.696175 2.457121 0.095131 0.248862 0.140825
0.218959 0.314440 0.612025 0.135107 1.165532 0.257336 0.120037 0.054679 5.306834 0.232523 0.299648 0.131932 0.481306 7.803902 0.089613 0.400547 0.245841 3.151815
2.547870 0.170887 0.083688 0.037967 1.959291 0.210332 0.245034 0.076701 0.119013 10.649107 1.702745 0.185202 1.898718 0.654683 0.296501 0.098369 2.188158 0.189510 0.249313

```

Where the amino acid frequencies are:

0.079066 0.055941 0.041977 0.053052 0.012937 0.040767 0.071586 0.057337 0.022355 0.062157 0.099081 0.064600 0.022951 0.042302 0.044040 0.061197 0.053287 0.012066 0.034155 0.069147

**Amino acid order** (matrix and frequencies)

A R N D C Q E G H I L K M F P S T W Y V

**Reference** Le S.Q., Gascuel O. Molecular Biology and Evolution. 2008 25(7):1307-20.
> 3 votes

---
Tags: phylogenetics, proteins, phylogeny
---


title: smallest length of reads in 3'quantseq?

text:
# Question
Title: smallest length of reads in 3'quantseq?
I am studying about RNA seq, and especially about 3'Quantseq(tagseq, 3primeseq). I wonder if there is a cutoff for the reads length. By this I mean that given that 3'quantseq targets the end of the genes and 3'UTR, then for example if reads are 75bp, then could be the case that only 3'UTR is mapped and not coding region? I am interested in Human Sapiens. Thank you in advance.

# Answer
Simple answer: Yes, that can very well happen depending on the average fragment size, the read length and the length of the UTR.
> 1 votes

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


# Question
Title: Problems replicating Kover paper results
I am trying to replicate the results for the creator of this repo: https://github.com/aldro61/kover2\_paper for genotype to phenotype machine learning interpretable (decision tree) prediction.

Specifically in the data section when I try to create the dataset as explained here: https://github.com/aldro61/kover2\_paper/tree/master/data

I make a directory called genomes\_dir in the data directory. Then I do:

```
python download_genomes.py 'mycobacterium tuberculosis/kanamycin/metadata.tsv' genomes_dir

```

at the end of the download it tells me that it saved file `genome_paths.tsv` the data directory, not genomes\_dir. So I run command:

```
kover dataset create from-tsv --genomic-data genome_paths.tsv ---phenotype-description "Kanamycin resistance" --phenotype-metadata 'mycobacterium tuberculosis/kanamycin/metadata.tsv  --output data_example.kover --progress

```

and it just hangs with a prompt:

```
(kover) nobu@gold3forever:~/Desktop/BioInformatics/ResearchAssistant/kover2_paper/data$ kover dataset create from-tsv --genomic-data genome_paths.tsv ---phenotype-description "Kanamycin resistance" --phenotype-metadata 'mycobacterium tuberculosis/kanamycin/metadata.tsv  --output data_example.kover --progress
> ^C

```

it does nothing after hours and hours so I control-C'd to end the process.

Does anyone have any idea what is happening?

I also had to install this to get the datasets from PATRIC: https://github.com/aldro61/patric\_tools

# Answer
The single quote at `'mycobacterium` is not matched so the command prompt is waiting indefinitely for the command to be completed. Maybe you want `'mycobacterium tuberculosis/kanamycin/metadata.tsv'` (note closing quote)
> 1 votes

---
Tags: machine-learning
---


title: Too slow issue of BioMart

text:
# Question
Title: Too slow issue of BioMart
I am using library(biomaRt) in R to annotate SNP rsid. It takes too much time for only 1000 SNPs....... I would like to know how I can run this code faster.

```
 coords<-data
 coords$chr_name<-coords$chromosome
 coords$start<-coords$position
 coords$end<-coords$position
 coords<-coords%>%dplyr::select(chr_name,start,end)
 mart=useMart(biomart="ENSEMBL_MART_SNP", host="grch37.ensembl.org",
 path="/biomart/martservice", dataset="hsapiens_snp")
 results=data.frame(matrix(ncol=4,nrow=0))
 for (i in seq (1,nrow(coords)))
{
ens=getBM(attributes = c('refsnp_id','allele','chrom_start','chrom_strand'),
filters = c('chr_name','start','end'), values = as.list(coords[i,]), mart = mart)
results=rbind(results,ens)}

```

I am using the following code, but still not working......

```
coords$newposition<-paste0(coords$chromosome,":",coords$position)
coords$chromosomal_region<-paste0(coords$newposition,":",coords$position)
coords<-coords$chromosomal_region
res<-getBM(attributes = c('refsnp_id','allele','chrom_start','chrom_strand'),
   filters = 'chromosomal_region',
   values = coords,
   mart = mart)

```

# Answer
Don't run getBM 1000 times. It can handle whole vectors, and get you all your results at once.
> 1 votes

---
Tags: r, biomart
---


# Question
Title: I want to display color key on my heatmap
```
library(gplots)

data=structure(c(-0.1152, 0.1548, 0.0037, 0, 0, 0, 0, 0, 0, 0.1369,
0.1914, 0, 0.0499, 0, 0, 0, 0, 0.1259, 0, 0.1854, -0.4412, 0,
-0.4425, 0, 0, -0.5141, 0, -0.3463, 0.1529, -0.1176, 0.4153,
-0.4412, 0.2169, 0.1894, 0.3288, 0.5549, -0.2549, 0.1915, -0.0759,
0.1255, 0.2269, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.4412, -0.1923,
0.1502, 0, 0, 0, -0.1978, 0.0828, 0, 0.2214, 0.192, 0, -0.2651,
0.192, -0.1412, -0.287, 0, -0.2299, 0, 0.1095, 0, 0.0166, 0,
0.1617, -0.16, 0, 0, 0.1255, 0, 0.2298, 0.153, 0.1627, 0, 0.293,
0.2476, -0.1565, 0.0899, 0, 0.2217, 0.2274, 0.1764, 0.1538, 0,
0.1747, 0.1857, -0.2309, 0.0827, -0.0609, 0.0754, -0.2653, 0.1461,
0.1611, -0.1201, 0, 0, 0.114, 0.237, 0.1738, 0, 0.2269, 0, 0.0068,
0, 0.1288, -0.1727, -0.1187, -0.1502, -0.1502, 0, -0.4104, -0.1305,
-0.1126, -0.2642, -0.1264, -0.2403, 0.1522, 0.0211, 0.1175, -0.0549,
-0.2571, -0.1305, -0.2143, -0.2153, -0.0566, -0.2068, 0.0406,
-0.087, 0.079, -0.0238), .Dim = c(10L, 14L), .Dimnames = list(
    c("AI", "AUDC", "AUDT", "AUT", "BD", "CB", "CD", "DPW", "MDD",
    "SZ"), c("BI", "BB", "CD", "CGC", "CTF", "HDT", "MTS", "NOI",
    "OFI", "SRV", "SAF", "FFI", "SUM", "WIN")))

data <- read.csv("C:/Users/Lenovo/Desktop/Correlation matrix VITD _Disorder_R/Z1.txt", header=TRUE, row.names=1,sep="")
data = as.matrix(data[![enter image description here][1]][1])
heatmap.2(data,
          trace ="none",
          density = "none",
          col = redblue(100),
         key=TRUE,
         keysize=1,
          cexRow =1,
          cexCol =1,
          margins = c(5,5),
          labCol = colnames(data),
          scale="row",
          hclust=function(x) hclust(x,method = "complete"))

  [1]: https://i.stack.imgur.com/jlwf4.png

```

# Answer
I'm not sure exactly what was going on with your data, but I read it in as a raw data.frame using `as.matrix(read.table(., row.names=1, header=T))`. The issue was that you were reading it in without a header, so the col.names were being read in as part of the matrix. But apart from that, the rest of your command was fine:

```
library(gplots)

data = structure(c(-0.1152, 0.1548, 0.0037, 0, 0, 0, 0, 0, 0, 0.1369,
0.1914, 0, 0.0499, 0, 0, 0, 0, 0.1259, 0, 0.1854, -0.4412, 0,
-0.4425, 0, 0, -0.5141, 0, -0.3463, 0.1529, -0.1176, 0.4153,
-0.4412, 0.2169, 0.1894, 0.3288, 0.5549, -0.2549, 0.1915, -0.0759,
0.1255, 0.2269, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.4412, -0.1923,
0.1502, 0, 0, 0, -0.1978, 0.0828, 0, 0.2214, 0.192, 0, -0.2651,
0.192, -0.1412, -0.287, 0, -0.2299, 0, 0.1095, 0, 0.0166, 0,
0.1617, -0.16, 0, 0, 0.1255, 0, 0.2298, 0.153, 0.1627, 0, 0.293,
0.2476, -0.1565, 0.0899, 0, 0.2217, 0.2274, 0.1764, 0.1538, 0,
0.1747, 0.1857, -0.2309, 0.0827, -0.0609, 0.0754, -0.2653, 0.1461,
0.1611, -0.1201, 0, 0, 0.114, 0.237, 0.1738, 0, 0.2269, 0, 0.0068,
0, 0.1288, -0.1727, -0.1187, -0.1502, -0.1502, 0, -0.4104, -0.1305,
-0.1126, -0.2642, -0.1264, -0.2403, 0.1522, 0.0211, 0.1175, -0.0549,
-0.2571, -0.1305, -0.2143, -0.2153, -0.0566, -0.2068, 0.0406,
-0.087, 0.079, -0.0238), .Dim = c(10L, 14L), .Dimnames = list(
    c("AI", "AUDC", "AUDT", "AUT", "BD", "CB", "CD", "DPW", "MDD",
    "SZ"), c("BI", "BB", "CD", "CGC", "CTF", "HDT", "MTS", "NOI",
    "OFI", "SRV", "SAF", "FFI", "SUM", "WIN")))

heatmap.2(data,
          trace ="none",
          density = "none",
          col = redblue(100),
         key=TRUE,
         keysize=1,
          cexRow =1,
          cexCol =1,
          margins = c(5,5),
          labCol = colnames(data),
          scale="row",
          hclust=function(x) hclust(x,method = "complete"))

```
> 1 votes

---
Tags: heatmap
---


title: error in random forest analysis

text:
# Question
Title: error in random forest analysis
I am now struggling to do random forest analysis, I will be thankful if you could help with code for random forest analysis.

I got samples from the root, soil, and leaf from two regions (bau & mau) and these samples belong to two seasons (Wet and Dry).

Now I am interested to do random forest analysis at genera or family level to identify the taxa which contribute the differences like in root samples based on region as well as season.

Here is my code, but I am getting the error.

```
    library(randomForest)
    library(knitr)

    #### RANDOM FOREST ANALYSIS #####
     #### Prepare data ####
     #Load OTU table
     OTU_table=t(read.table("asv.table.txt", row.names=1,sep="\t", header=T, blank.lines.skip=F, check.names=F))
     table(apply(OTU_table,1,sum)) #verify rarefaction

     #Load metadata
     Meta=read.table("metadata.txt", header=T, row.names=1, stringsAsFactors=F, na.strings="NA",check.names=FALSE)
     Meta$sampleid=rownames(Meta)

     #Load taxonomy
     taxo=read.table("taxonomy.txt", row.names=1, sep="\t", header=F ,stringsAsFactors=F,quote="")
     rownames(taxo)=paste("a.",row.names(taxo),sep="")

     #### Run models ####
     #1. Root only
     # 1.1. both region
     # 1.2. bau only
     # 1.3. mau only
     #2. Soil only
    # 2.1. both region
     # 2.2. bau only
     # 2.3. mau only
     #Params RF
     NTREE=1000 # Number of Trees
     NbVar=1000 # Number of variables tested at each split

     #### Root ONLY 1-3 ####

 # 1. BOTH Region

 #Subset of data
 RootSamples=as.character(Meta[Meta$Compartment=="Root","sampleid"])
 Root_OTU_table=OTU_table[RootSamples,]

 #Model with microbiome based on Season, region

whole_root_pred=data.frame(Season=Meta[RootSamples,"Season"],Region=Meta[RootSamples,"Region"],a=Root_OTU_table)

head(whole_root_pred)
     Season Region a.d2ec9f3b77975c0f457e4b7413b217ff
     a.3147790f0d5a78316fb9dd64f53b9473 a.97aecc1f35cc1f50db507ad71dd22367
     a.bfad6370d28182cc6304844e9bec7fb6 a.5fa2a987221a1d9ca416148570c18086

    **RF_model_Root_all=randomForest(y=?,sampsize=c(143,143),strata=?,x=whole_Root_pred,importance = T,proximity = T,ntree =
    NTREE,mtry = NbVar)**

print(RF_model_Root_all)
     #plot summary using the 5% most important OTUs ERROR ON LAST LINE
     imp=data.frame(importance(RF_model_Root_all))
     imp$genus=as.character(taxo[rownames(imp),"Genus"])
     Best=imp[imp$MeanDecreaseAccuracy>quantile(x = imp$MeanDecreaseAccuracy,.95),]
     bymedian <- with(Best, reorder(genus, -MeanDecreaseAccuracy, median))

     pdf(width = 20,height = 10,file=paste(pathforplots,"Variable_Importance_Root_BothRegion_raref.pdf",sep=""))
     par(mar=c(15,5,1,1))
     boxplot(Best$MeanDecreaseAccuracy ~ bymedian, data = Best,
     xlab = "", ylab = "Variable Importance",
     main = paste("Root in Both Countries; Error Rate=",round(RF_model_Feces_all$err.rate[NTREE,"OOB"],3),sep=""), varwidth = TRUE,
     col = "lightgray",las=2)

     dev.off()

```

Many thanks

# Answer
> 2 votes
This question will be a little hard to answer without more information.

For example, we will need to see your dataset (`whole_root_pred`), to decide why `Stunting_Root` is NULL.

1. You might need to initialize `Stunting_Root` as a variable. It is currently not clear if it is e.g. a column of your dataframe, or just uninitialized. Uninitialized variables are NULL, which would explain your problem. randomForest might not know to look for `strata` inside your dataframe, for example. Is it in your dataframe?
2. Also, I might be missing something, but why are you passing `?` as a response? I'm not an expert but I believe that is an illegal character in R (I'm pretty sure?).

---
Tags: r, bioconductor, statistics, microbial-genomics
---


# Question
Title: panic: division of zero by zero or infinity by infinity error when running GODON positive selection analysis
I am using `godon` to find genes under positive selection with the **branch-site test**. However, after running for a little while I get the error *"panic: division of zero by zero or infinity by infinity"*.

Below is the full output with the error report. I can supply the phylogenetic tree and fasta file I am using if that would be helpful. Thank you.

> iter=33 lnL=-73104.072
> iter=66 lnL=-72956.855
> iter=99 lnL=-72912.956
> iter=131 lnL=-72902.900
> iter=162 lnL=-72900.147
> iter=194 lnL=-72898.562
> iter=226 lnL=-72897.970
> iter=258 lnL=-72897.651
> iter=289 lnL=-72897.435
> iter=321 lnL=-72897.284
> iter=353 lnL=-72897.199
> iter=385 lnL=-72897.153
> iter=416 lnL=-72897.108
> iter=448 lnL=-72897.076
> iter=17 lnL=-72890.551
> panic: division of zero by zero or infinity by infinity goroutine 1 \[running\]: math/big.(\*Float).Quo(0xc4273a8f40, 0xc4273a8f68, 0xc4273a8f18, 0xc42e88c118) /usr/local/go/src/math/big/float.go:1639 +0x15c bitbucket.org/Davydov/godon/cmodel.(\*BranchSite).BEBPosterior(0xc422164120, 0xc4259bb900, 0x91, 0x91) /go/src/bitbucket.org/Davydov/godon/cmodel/bs.go:416 +0x6e9 bitbucket.org/Davydov/godon/cmodel.(\*BranchSite).Final(0xc422164120, 0x101) /go/src/bitbucket.org/Davydov/godon/cmodel/bs.go:455 +0x112 main.computeFinal(0x1f263a0, 0xc422164120, 0xc4201522a0, 0x101, 0x6a80, 0x0) /go/src/bitbucket.org/Davydov/godon/godon/optimize.go:150 +0x93 main.performSingleTest(0xc420162320, 0xc428266000, 0x9399, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, ...) /go/src/bitbucket.org/Davydov/godon/godon/hyptest.go:335 +0x1b52 main.hypTest(0x1ee03a7, 0x1ee03a7, 0x4, 0xc4201aa501, 0x1, 0x1) /go/src/bitbucket.org/Davydov/godon/godon/hyptest.go:88 +0x71a main.main() /go/src/bitbucket.org/Davydov/godon/godon/godon.go:316 +0xefa

# Answer
> 1 votes
Ah wait you mentioned a branch sites test, right that makes sense. You've got two or more identical sequences in the analysis. Ensure there are no identical sequences in your alignments and ... PRESTO! It will work :-)

Below are my original thoughts before I re-read your question. I've kept the notes because it helps explain the theory a bit. There is no need to supply the tree, it should be resolved.

---

The tree likelihood looks to be decreasing nicely, so it can see both the tree and the codons in the alignment.

At a guess the error relates to two or more indentical nucleotide sequences in your alignment. Thus zero non-synonymous sites (amino acid change) and zero synonymous (silent) sites. This error would only occur if you were calculating selection for each 'node' in the tree (internal branch). Alternatively, possibly if the number of synonymous sites was zero, regardless of the non-synonymous mutations this error could appear, but this doesn't seem to fit the error output.

I don't know "Godon", but recommend `PAML` for this analysis. Admittedly it is complicated to deploy, but it is a very robust algorithm and hugely respected.

---
Tags: phylogenetics, phylogeny, positive-selection, codon, godon
---


title: BWA doesn't recognize HiC options

text:
# Question
Title: BWA doesn't recognize HiC options
I'm trying to run bwa on HiC data. This is the command I'm giving:

```
bwa mem -5SP -T0 -t4 -R '@RG\tID:\tSM:\tLB:\tPL:ILLUMINA\tPU:none' REFERENCE.fa PAIR1.fq PAIR2.fq

```

and here's the error message:

```
mem: illegal option -- 5

```

Supposedly, the 5SP option is for HiC data (cite), but my BWA isn't recognizing it. Anybody have any ideas? My initial thought is to update BWA to the newest version, but I don't want to mess with our replicatibility since I am new to the lab and don't want to be updating software when I don't need to. The server is a Mac Pro and bwa version is 0.7.12-r1039. Thanks.

# Answer
> 3 votes
This option has been added in 0.7.16a, see https://github.com/lh3/bwa/releases, so consider updating.

# Answer
> 3 votes
Your link with those command line options is 1 year old. Looks like that version of bwa is 6 years old. You probably should update.

---
Tags: bwa, hi-c
---


# Question
Title: Looking for genes that are expressed the same in all the samples
I have 5 samples (from a RNA seq experiment) from patients but I have no control, so since we can't look for the differentially expressed genes we want to look for the genes that have similar expression in all the sample. We want to do that to see if some genes that can be involved are among those expressed in all the samples. I have never done that, I do not know how to do it and I cant find any info about it, so any ideas are welcome.

Thanks!

# Answer
I cannot think of a way without performing differential expression analysis. Since you don't have controls, what you could do is assigning some of the samples as controls.

Let's say you have samples A, B, C, D, E and F. You can perform differential expression analysis A, B, C (i.e. treatment) versus D, E (control). If you perform this for all possible combinations (that would be 5 choose 2 which is equal to 5 choose 3 if you decide to go for 3 vs 2 splits) and filter for genes which do not show up in any DE list you end up with the genes that are expressed similarly in al of these samples.
> 1 votes

---
Tags: gene-expression
---


# Question
Title: install bowtie2 from sources cannot find -ltbb
I am trying to install bowtie2 aligner from sources without the root access. Bowtie2 needs tbb package to be installed and it is recommended to install oneTBB.

```
# install one TBB
mkdir -p SHOME/soft
cd $HOME/soft
git clone https://github.com/oneapi-src/oneTBB.git
cd oneTBB
make

```

Then I add C\_INCLUDE\_PATH to .bashrc:

```
echo "export C_INCLUDE_PATH=$HOME/soft/oneTBB/include" >> ~./bashrc

```

and continue to install bowtie2:

```
# install bowtie2
mkdir -p $HOME/soft/bowtie/source
cd $HOME/soft/bowtie/source
wget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.4.1/bowtie2-2.4.1-source.zip
unzip bowtie2-2.4.1-source.zip -d ..
cd ../bowtie2-2.4.1
make

```

TBB library is recognized but in the end I get this error:

```
/usr/bin/ld: cannot find -ltbb
/usr/bin/ld: cannot find -ltbbmalloc_proxy
collect2: error: ld returned 1 exit status

```

ltbb and ltbballoc\_proxy are not in TBB but libtbb and libtbballoc\_proxy instead. I also tryed to add variables to .bashrc:

```
TBB_LIB_PATH=$HOME/soft/oneTBB
TBB_INCLUDE_PATH=$HOME/soft/oneTBB/include
TBB_RELEASE_DIR=$HOME/soft/oneTBB/build/home/suvar/apps/oneTBB/build/linux_intel64_gcc_cc9.3.0_libc2.31_kernel5.4.0_release

```

but seems it does not help.

I need help to complete the installation of bowtie2 on our server.

# Answer
You need to give the path to the source files for the compiler *and* the full path to the compiled libraries (all parent directories included) for the linker. You also need to use either CPATH or CPLUS\_INCLUDE\_PATH instead of C\_INCLUDE\_PATH to point to the oneTBB source, since bowtie2 is C++ and C\_INCLUDE\_PATH will only have an effect for C (see Environment Variables Affecting GCC in the GCC documentation for more info.)

For example like this, inside bowtie2-2.4.1, without making any lasting changes to the shell's setup in this case:

```
$ ONETBB=$HOME/soft/oneTBB
$ export LIBRARY_PATH=$ONETBB/build/linux_intel64_gcc_cc5.4.0_libc2.23_kernel4.4.0_release
$ export LD_LIBRARY_PATH=$LIBRARY_PATH
$ export CPATH=$ONETBB/include
$ make

```

After it's built you'll still need the library path to oneTBB defined (via LD\_LIBRARY\_PATH) so bowtie2 can find it at run-time.

**But if you don't have a specific reason to compile from source I'd recommend avoiding it.** You could go with what h.mon suggested and use an existing binary, or use a package manger to install it like zorbax suggested. Conda has a package for bowtie2 in the bioconda channel, up through and beyond version 2.4.1, and you can use that as a package manager without root access.
> 1 votes

---
Tags: software-installation, bowtie2
---


# Question
Title: Import a tab-separated file with differing numbers of elements in each row; prokka output
I am using prokka to annotate a bacterial genome:

`prokka ecoli.fa`

Prokka is outputting a tab-separated file (called PROKKA\_12142017.tsv) with differing numbers of elements in each row:

```
locus_tag    ftype    gene    EC_number    product
BOFHCHGE_00001    CDS    hypothetical protein
BOFHCHGE_00002    CDS    rdgB_1    3.6.1.9    dITP/XTP pyrophosphatase
BOFHCHGE_00003    CDS    hypothetical protein
BOFHCHGE_00004    CDS    hypothetical protein
BOFHCHGE_00005    CDS    hypothetical protein

```

Does anybody know how I can get round this to import into a dataframe in R?

Is there a way to ask prokka to add something in the gene column?

# Answer
Author here - this was a bug which was fixed in August:

See the issue in Github.

You will need to use the latest Github HEAD version until the next release.
> 5 votes

# Answer
You don't really need to get around anything, your example can be loaded correctly with `read.delim()`. The rows with missing values are filled in with blanks. You may, however, prefer the `readr` package, which handles this a bit more elegantly (it'll tell you which values were missing and not do the annoying string-\>factor conversion).
> 6 votes

# Answer
As Devon said, `read.delim` can deal with this perfectly well as long as the file is properly formatted (so missing fields still end with the field separator, like `\t\t`). For example, with this input:

```
Field1  Field2  Field3
Row1field1  Row1field2  Row1field3
Row2field1      Row2field3

```

Or, to show the fields more clearly:

```
Field1\tField2\tField3
Row1field1\tRow1field2\tRow1field3
Row2field1\t\tRow2field3

```

R can import:

```
> df<-read.delim("file",header=TRUE, sep="\t")
> df
      Field1     Field2     Field3
1 Row1field1 Row1field2 Row1field3
2 Row2field1            Row2field3

```

Alternatively, if you really want to add something to the missing columns, you can use:

```
$ awk -F"\t" '($3==""){$3="none" }1;' file
locus_tag   ftype   gene    EC_number   product
BOFHCHGE_00001  CDS hypothetical    protein
BOFHCHGE_00002  CDS rdgB_1  3.6.1.9 dITP/XTP    pyrophosphatase
BOFHCHGE_00003  CDS hypothetical    protein
BOFHCHGE_00004  CDS hypothetical    protein
BOFHCHGE_00005  CDS hypothetical    protein

```

Or, to add `none` to any empty field, not just the 3rd one:

```
awk -F"\t" '{for(i=1;i<=NF;i++){if($i==""){$i="none"; }}}1;' file > newfile

```
> 3 votes

# Answer
Since the missing values don't have field separators, the easiest option is to rearrange the columns so that the ones with missing values are at the end. They will then be the ones that are filled with missing values when you import them into R. You can use use `fill=TRUE` with read.table() or use readr, which will complain about missing columns, but fill them in with missing values automatically:

```
awk -F'\t' -v OFS='\t' '{print($1,$2,$5,$3,$4)}' file.txt |perl -pe 's/\t+/\t/g' >newfile.txt

```

With readr:

```
read_tsv("newfile.txt")

# A tibble: 5 x 5
       locus_tag ftype                  product   gene EC_number
           <chr> <chr>                    <chr>  <chr>     <chr>
1 BOFHCHGE_00001   CDS     hypothetical protein   <NA>      <NA>
2 BOFHCHGE_00002   CDS dITP/XTP pyrophosphatase rdgB_1   3.6.1.9
3 BOFHCHGE_00003   CDS     hypothetical protein   <NA>      <NA>
4 BOFHCHGE_00004   CDS     hypothetical protein   <NA>      <NA>
5 BOFHCHGE_00005   CDS     hypothetical protein   <NA>      <NA>

```

The perl regex is necessary because awk adds tabs after the missing values, which defeats the point of reordering the columns
> 1 votes

# Answer
all. I know this question was asked 2 years ago, but I wrote a script in R, with a little bash, to parse the tsv file in the old version of prokka. In the new version of prokka this problem is fixed.

First use bash to get the rows where ftype is CDS.

```
grep "CDS" all_prokka_annotations.tsv > all_prokka_annotations_CDS.tsv

```

Next we can parse the TSV with the CDS lines in R.

prokka_annotations = readLines("all_prokka_annotations_CDS.tsv")
prokka_annotations_split = strsplit(prokka_annotations,split="\t")
# first get the protein product name.
prokka_annotations_product = sapply(prokka_annotations_split, function(x) x[length(x)])
# now the I have all protein product names remove them
prokka_annotations_split = sapply(prokka_annotations_split, function(x) x[-length(x)])
#next get locus tag and ftype
locus_tag_cds = sapply(prokka_annotations_split, function(x) x[c(1,2)])
locus_tag_cds = t(locus_tag_cds)
# now remove locus tag and ftype
prokka_annotations_split = sapply(prokka_annotations_split, function(x) x[-c(1,2)])
# next get EC number
ec_number = lapply(prokka_annotations_split, function(x) x[grep(".",x,fixed=T)])
ec_number = as.character(ec_number)
ec_number[ec_number=="character(0)"] <- NA
## now get the gene names if available. to do this remove the EC numbers. whats left are gene names
prokka_gene_names = lapply(prokka_annotations_split, function(x) x[!grepl(".",x,fixed=T)])
prokka_gene_names = as.character(prokka_gene_names)
prokka_gene_names[prokka_gene_names=="character(0)"] <- NA
prokka_parsed_tsv_correct = cbind(locus_tag_cds,product=prokka_annotations_product,gene=prokka_gene_names,EC=ec_number)

```

> 1 votes

---
Tags: r, sequence-annotation
---


# Question
Title: Trimmomatic trims most of the reverse strand
I have a paired-end data. It is from small part of human genome not a WGS or WES. I use BWA for alignment and do variant calling. Since I had some false positives I wanted to do trimming.

I used trimmomatic with command:

```
trimmomatic PE -phred33 -threads 2 read1.fastq.gz read2.fastq.gz read1_trimmed.fastq read2_trimmed.fastq LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15`

```

I do not know why but I have this output:

```
Input Read Pairs: 1834513
Both Surviving: 1822034 (99.32%)
Forward Only Surviving: 5740 (0.31%)
Reverse Only Surviving: 417 (0.02%)
Dropped: 6322 (0.34%)
TrimmomaticPE: Completed successfully

```

I do not have any adapters since the are removed in `bcl2fastq` step. In FASTQ quality control I have seen that quality of reverse and forward reads are similar in quality but trimmomatic trimmed many reads in reverse file. I cannot do alignment with these files, BWA gives a warning that says:

```
[bseq_read] the 2nd file has fewer sequences.
[process] read 155012 sequences (16693324 bp)...
[bseq_read] the 2nd file has fewer sequences.

```

Is there anyone who faced a problem like this before?

Thank you in advance

# Answer
> 1 votes
Maybe you need to remove the unpaired reads:

```
trimmomatic PE -phred33 -threads 2 \
    read1.fastq.gz read2.fastq.gz \
    read1_trimmed.fastq.gz read1_unique.fastq.gz \
    read2_trimmed.fastq.gz read2_unique.fastq.gz \
    LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15

```

Then you can try with bwa.

---
Tags: ngs, quality-control, trimming
---


title: WGBS vs ONT, which one should I trust?

text:
# Question
Title: WGBS vs ONT, which one should I trust?
First of all let me apologize if this is not the place for this question. While it can be quite broad, I will try to make it as specific as I can.

**Context**

I am testing the correlation between Whole Genome Bisulfite Sequencing and Oxford Nanopore methylation calls for one individual (mammal, blood tissue) with replicates at all CpG sites. In the (very) best cases I find a close to 0.9 correlation. Yet, for the part that does not correlate I would like to know which technology is closer to the "biological truth".

**My strategy**

My plan would be to first find CpG sites/short CpG regions showing differential methylation states between WGBS and ONT. Then I would validate the methylation states of these sites using a third party technology such as pyromark.

**What I am asking to the SE community**

* First, has anyone done this before ? I would love to hear about your experience
* Do you see some critical aspects I should be watching out for ?
* Last if you have comments about the strategy I propose here, you are welcome

# Answer
Methylation can be cell-specific, which makes it difficult to evaluate accuracy on a bulk-cell level (even within the same tissue). How can you tell that the differences you're seeing are due to platform differences, or due to biological variation?

I find that adding more haystacks doesn't help much in working out the truth of a dataset. If you want to investigate biological truth, it would be better to create a biological system with a known methylation (or demethylation) pattern and test that.
> 2 votes

# Answer
I would assess directionality and accuracy of prediction by 1) WGBS predicting ONT and then 2) ONT predicting WGBS.

Firstly, I would use deep learning (or machine learning) and train WGBS against ONT, parameterise and then test. Then conversely train ONT against WGBS, parameterise and then test. The approach with highest accuracy of prediction (using the 'accuracy' index) would be the approach assumed to be "closer to the biological truth". If both calculations produced comparable accuracies I would conclude what @gringer has stated that natural variation / natural heterogeneity is the predominant signal in the signal. This conclusion depends on the ability of deep learning to map the biological process.

It is an approach highly applicable to a deep learning estimate regardless and may not replace a truely controlled test, this will depend on what controls you have run, but could provide a valuable and easily obtainable insight in its absence.

If you did run a control sample, then you are truely in, because this would provide the training set and the WBGS and ONT provide the test data.

---

You asked me how it works It doesn't replace good controls, but it is 'trendy workaround'.

Two steps

1. In deep learning (artificial neural networks, ANN, in this case you would use RNN such as LSTM) or machine learning (ML) e.g. random forest - random forests are MUCH easier to do, you supply a series of known values, thus CpG sites and their WGBS scores and use ANN (RNN approach) or ML (e.g. random forest). Both methods will look for pattens between the observed DNA data and the CpG occurence (training set). Usually this is 60% of your data set is used in training. A separate 20% of the data set is often used for parameterisation and a final 20% of the data set is used for testing. Thus how well can a ANN or ML predict a CpG site based on WGBS data. Likewise you would do the same for ONT. There is a caveat here which I discuss below.

The idea is if you give it a bit of DNA it will assign whether it is CpG (see caveat below).

2. However, instead of testing WGBS to WGBS you also test it WGBS to ONT and ONT to WGBS and the method with the highest 'accuracy' (positive hits, against true hits) is 'better method'. This is because you already know which is ONT positive and ask a WGBS trained algorithm whether it thinks its positive. The relationship between observed CpG and 'true (experimental) CpG is your accuracy score.

The caveat is that you need to supply the training set with sites that are known to be negative under WGBS or ONT. You need to think the best strategy for doing this. If some of those negatives were positve in the other method that would affect the outcome alot.

There is an issue about vectorisation this is often done using k-mers, and is a headache. Basically you have to make a bit of DNA into numbers and that in my opinion is the difficult part. In your case I ain't so sure its so difficult. If you can get that bit correct and biologically meaningful you are in business. My stuff is phylogeny based and I understand the relationship between biology (mutation) and numbers. In your case that needs thought, however someone will have thought about it and solved it.

As a personal statement I would do this via random forests, because ANN is difficult.

One such vectorisation method is here by Li et al (2017), however as they have used it, the method has flaws due to equal weighting of mutation frequencies. However, these flaws do not apply to your work - it would actually be a good method, they singly apply to tree building.

---

The final thing is that the more experimental data you can feed it the better. You don't have to complicated about it, just create a new column in the input section, e.g. along side teh vectorised DNA and the algorithm will figure out whether it helps. You can't feed it enough data, because it is doesn't help it wouldn't use it.
> 1 votes

---
Tags: methylation, epigenetics
---


title: dbnSNP frequency anomalies

text:
# Question
Title: dbnSNP frequency anomalies
Sometimes dbSNP reports very different allele frequencies for different large-scale genome projects e.g. between 1000 Genomes and GnomAD

rs11822440 1000Genomes A=0.4629 C=0.5371 GnomAD A=0.99997 T=0.00003

rs28531832 1000Genomes T=0.4379 C=0.5621 GnomAD T=0.99993 A=0.00007

rs78931618 1000Genomes A=0.7334 C=0.2666 GnomAD A=0.9999 T=0.0001

rs75716913 1000Genomes T=0.9006 C=0.0994 GnomAD T=0.99149 A=0.00851

Is this an issue with one or other of these projects data processing? Is something else going on here biologically or technically that might cause this?

# Answer
> 1 votes
I'll stick my neck out and simply state that allopatry alone would account for the pattern you are observing. Thus populations can exihibit variation in allelic profiles, particularly where Mendellian inheritance is operating (i.e. higher eukaryotes), singly due to geographic isolation.

This represents the Meyers view of population genetics and would account for variation within neutral theory, e.g. silent mutations. Obviously I'm guessing the populations were selected by geographic foci in generating the allelic frequencies.

---
Tags: snp, 1000genomes
---


# Question
Title: How to subset a BAM by a list of QNAMEs?
I have a text file 'qnames.txt' with QNAMEs in the following format:

```
EXAMPLE:QNAME1
EXAMPLE:QNAME2
EXAMPLE:QNAME3
EXAMPLE:QNAME4
EXAMPLE:QNAME5

```

I would like to subset my BAM `file.bam` via all of these QNAMEs into a new SAM.

Naturally, I can do this individually, e.g.

```
samtools view file.bam | grep 'EXAMPLE:QNAME1' > subset.bam

```

But I'm unsure how to do this for a list of QNAMES:

1. How do I write a for loop that will do all of these queries, outputting the correct SAM needed?

I could write a `for` loop that creates *n* SAM files, and then `cat`s them…
2. Is there a way to specifically `grep` by QNAME? The above may `grep` reads that may not be associated with the correct QNAME?
3. How do I keep the BAM header?

# Answer
(1) It won't be super fast but you can provide grep with a file of QNAMES.

```
samtools view file.bam | grep -f 'qnames.txt > subset.sam

```

where qnames.txt has

```
EXAMPLE:QNAME1
EXAMPLE:QNAME2
EXAMPLE:QNAME3
EXAMPLE:QNAME4
EXAMPLE:QNAME5

```

(2) This would be a little more complicated but can you give an example where the grep might be have the correct QNAME?

(3) To keep the BAM header I would use `samtools -H file.bam > header.txt` to get the header and then cat the header file with the grep'ed file
> 4 votes

# Answer
use picard FilterSamReads http://broadinstitute.github.io/picard/command-line-overview.html#FilterSamReads

> READ\_LIST\_FILE (File) Read List File containing reads that will be included or excluded from the OUTPUT SAM or BAM file. Default value: null.
> 9 votes

# Answer
**UPDATE 2021/06/28:** since version 1.12, samtools now accepts option `-N`, which takes a file containing read names of interest.

Using samtools 1.12 or greater:

```
samtools view -N qnames_list.txt -o filtered_output.bam input.bam

```

**OLD ANSWER:**

When it comes to filter by a list, this is my favourite (much faster than grep):

```
# given a bam file "aln.bam" and a list of read names "reads.txt":
samtools view -h aln.bam \
| awk 'FNR==NR {reads[$1];next} /^@/||($1 in reads)' reads.txt - \
| samtools view -b - > filtered_aln.bam

```

Basically, the awk command first parses the reads.txt file, creates an array of QNAMEs, then parses the alignments and print if the QNAME is present in the array (or if the line belongs to the header).
> 7 votes

# Answer
Below is a small bit of python code showing one naive but manageable way of doing this (N.B., I expect `grep` to be faster, though getting it to output the header will be annoying):

```
import pysam
qnames = set(...)  # read names go here
bam = pysam.AlignmentFile("some input file.bam")
obam = pysam.AlignmentFile("some output file.bam", "w", template=bam)

for b in bam.fetch(until_eof=True):
    if b.query_name in qnames:
        obam.write(b)
bam.close()
obam.close()

```

Now that'll work, but if you have a lot of names to look through then it'll end up being quite slow (as will `grep -f ...`). If you need to select a large list of reads by name, then a more efficient strategy is:

1. Name sort the BAM file (`samtools sort -n` or picard), which can be done with multiple threads.
2. Sort the list of read names to match (this will be more complicated than one would naively think, since samtools and picard will name sort differently, so be sure to put some thought into this).
3. Perform the same sort of iteration as above, but only compare to the top read name in your list from step 2 (removing that element from the stack after a match).
> 4 votes

# Answer
Performing a grep on $n$ alignments to $m$ qnames would give you $O(mn)$ operations. If you decide to sort your BAM and qnames first, you can reduce it to $O(\min\\{m, n\\})$ by popping reads and qnames off their stacks. You also don't read the entirety of the BAM at once this way, limiting memory consumption. I'm pretty sure Picard does something like this, as well.

Here's an example in Python that I've written for this very purpose:

```
import pysam

def filter_qname(bamfile, idfile, outfile):
    '''
    Filter reads in a SAM/BAM file by their query names

    Parameters
    ----------
    bamfile : str
        Sorted BAM file path
    idfile : str
        Text file path containing qnames to keep
    outfile : str
        Output file to write to
    '''
    # read in name-sorted BAM file
    bam = pysam.AlignmentFile(bamfile, "rb")
    if outfile.endswith('.sam'):
        output = pysam.AlignmentFile(outputfile, "w", template=bam)
    elif outfile.endswith('.bam'):
        output = pysam.AlignmentFile(outputfile, "wb", template=bam)
    else:
        raise ValueError("Unknown output file format for `outfile`: {}".format(outfile))
    #  read in IDs to be removed (list(set(...)) to only keep unique IDs)
    ids = list(set([l.rstrip() for l in open(idfile, 'r').readlines()]))
    # sort IDs to match BAM for efficient processing (destructive function)
    ids.sort(key=str.lower, reverse=True)
    n_ids = len(ids)
    # variable for ensuring BAM is sorted by query name
    last_q = None
    reads = bam.fetch(until_eof=True)
    read = next(reads)
    while True:
        # if read name is greater than current top of stack
        if read.query_name > ids[-1]:
            # if this is the last ID
            if n_ids == 1:
                # write remaining reads to new file and break out of while loop
                output.write(read)
                for r in bam:
                    output.write(read)
                break
            # otherwise pop that ID, try again
            else:
                ids.pop()
                n_ids -= 1
        # skip reads that match the top of the stack
        elif read.query_name == ids[-1]:
            read = next(reads)
        # if read name is less that top of stack, write it and move on
        else:
            output.write(read)
            read = next(reads)
    output.close()

```

I've implemented it in my own miscellaneous bioinformatics tools package here if you want to see it in action
> 3 votes

# Answer
Most of the solutions posted thus far are either slow, or can yield spurious results in certain cases. (I haven't tested Picard, which I assume will work as intended. But personally I tend to recoil when I have to set up a JVM!)

Let's approach the problem from a different angle:

Good hash functions have ~O(1) complexity for a lookup. If we can create a hash table of our QNAMEs, then searching a BAM file should be O(n) in the number of records. The Python developers are very proud of the hash function in cpython, so let's use that:

```
#!/usr/bin/env python3

import sys

with open(sys.argv[1], 'r') as indexfile:
    ids = set(l.rstrip('\r\n') for l in indexfile)

for line in sys.stdin:
    qname, _ = line.split('\t', 1)
    if qname in ids:
        sys.stdout.write(line)

```

With a test set consisting of a ~3GB BAM file and a `qnames.txt` with ~65K entries, running `samtools view input.bam | ./idfilter.py qnames.txt` takes roughly **6.5 seconds** on my test server. By comparison, piping the SAM to `fgrep -f qnames.txt` (using GNU grep) takes about 8.5 seconds (but may yield spurious results).

(Why might `grep` yield spurious results? If you have a QNAME `foo` in your `qnames.txt`, but your BAM file also contains `fooX` and `fooY`, they will all be matched by `fgrep`. The solution is to anchor each of your search patterns between the beginning of the line and a tab, e.g. `^foo\t`, but then you have to use standard `grep` which will have to construct an NFA for each pattern, and not `fgrep` which implements Aho-Corasick, and you'll be orders of magnitude slower.)

I've rewritten my little Python script in Rust using `std::collections::HashMap` for the hash function, and using the rust\_htslib crate to read and write straight from BAM, and this is more than twice as fast as the Python, even when reading and writing BAM. However, my Rust isn't really fit for public consumption…
> 2 votes

# Answer
You can use SAMsift:

```
samsift -i file.bam -0 'q=open("qnames.txt").read().splitlines()' -f 'QNAME in q'

```

`-i file.bam` specifies the input BAM file. `-0 'q=open("qnames.txt").read().splitlines()'` loads your list of QNAMEs. `-f 'QNAME in q'` specifies the filter for alignments.

With SAMsift one can use any Python expression for filtering. The downside is that it's slower than the other tools.
> 1 votes

---
Tags: sam, samtools
---


# Question
Title: Performing PCA for the samples and for the genes
I have 10 samples from a RNAseq experiment (5 control, 5 disease), I have performed a cluster analysis for the samples and for the genes (4000 genes aprox) to see how they cluster (to see which samples are similar and which genes have similar expressions). I was wondering if in terms of statistics it makes any sense to perform a PCA of the samples instead of the genes?

thank you in advance!

# Answer
First of all, I am not a statistician but have been performing RNA-seq analysis for a while, here is my take:

* I think performing PCA on the samples makes sense in terms of mathematics. What I mean is you can compute eigenvectors and subsequently principal components either on your genes or on your samples.
* I don't think performing PCA on the samples makes sense in terms of statistics in RNA-seq context (and probably in any other context): Your genes (features) describe your samples and not the other way around. So the computed PCs should describe/explain your samples and not the other way around. You would not perform PCA on your samples (individual cells in scRNA-seq) even in the case of scRNA-seq, where number of samples sometimes exceeds the number of genes, there PCA is performed on genes and the resulting few PCs are used for downstream analysis, usually clustering and visualization, both used to infer about samples (cells).
> 2 votes

# Answer
If you use, say, R's prcomp, you can use `my_pca_object$x` to get PC coordinates for the samples, and `my_pca_object$rotation` to get how much each gene contributes to each PCA. (You might need to transpose your transformed count file to get it to work right)
> 1 votes

# Answer
It is quite possible to make a PCA plot with the samples. Follow this tutorial : https://hbctraining.github.io/DGE\_workshop/lessons/03\_DGE\_QC\_analysis.html
> 1 votes

---
Tags: pca
---


title: .cif .pdb conversion with python

text:
# Question
Title: .cif .pdb conversion with python
can you please give me an advice how to convert .cif files into .pdb preferably using python?

Thanks in advance! Best, Balint Biro

# Answer
As @Michael said PyMOL can do this. PyMOL is not only an app, but also a python package —installed via conda not via the regular download (you can have both).

```
import pymol2
with pymol2.PyMOL() as pymol:
     pymol.cmd.load(infile,'myprotein')
     pymol.cmd.save(infile.replace('.cif', '.pdb'), selection='myprotein')

```

This runs the parallelisable instance, so you can have that context manager block running on how many threads or processes you fancy.
> 2 votes

---
Tags: python, pdb
---


title: Power law distribution alpha values

text:
# Question
Title: Power law distribution alpha values
I did a Pearson co-expression analysis for generating networks for my tissue-specific (chondrocytes) RNA-seq data. I used R package poweRlaw to check for the power-law distribution. We got the following values:

```
 Xmin = 14
 goodness of fit = 0.035
 α = 4.76

```
1. My question is what does alpha (α) of more than 4 indicates?
2. Is it is okay to have such values?

Thank you

# Answer
It means the distribution is heavily skewed and has a 'long-tail', thus the variance is much greater than the mean. In other words most values are around zero, with a small number of values having extreme values.

\[

You need to check the goodness of fit because that is saying the deviation from a power-laws distribution is just significant. So it ain't quite power-laws, but the general gist of a 'long-tailed distribution' is correct.

Its a good result and a nice analysis, but you are not quite there at present (needs to be \>0.05). Check the negative binomial distribution.

---

I've seen you data, it looks very good and the power series is not a bad explanation of the data. There are two possibilities for the goodness of fit being 0.035

1. The goodness of fit is low because of 'outliers' at the beginning of the distribution
2. An alternative distribution, such as a negative binomial distribution fits better
3. Check it is 0.035 not 0.35 because that changes everything.

The best way to do this is to plot the power series expected distribution onto the graph and look where the fit 'fails'. If it is not fitting at the beginning of the distribution thats is about outliers and easily solvable. If it is not really capturing the shape then its probably a slightly different distribution.

If its point 1 then just remove the very first point in the graph and try again and check the p-value. Or take the mean of the clump of three points next to the first point and try again (or just remove those three point). Don't do any more than that and the technical word is "outlier" analysis. In your write-up you must explain the issue around any 'outliers' and that is perfectly acceptable. If it is a power-series distrupted by an outlier(s) there is formal analysis to demonstrate an 'outlier' using outlier analysis if any one ever asks you.

If its point 2. then use e.g. negative bionomial or the gamma distribution (alpha \<1) and check the goodness of fit.

You want a goodness of fit \> 0.05 and preferably \>0.1 to be sure its the right distribution. There are formal statistics to do this, but that is the main one.

My guess is there are outliers at the beginning because your sample size is not large enough to identify the variation correctly, hence the fit is low. In 'long-tails' its easy to have a 'low sample size' because of the way they work. Thus in this case mapping a distribution is a perfect way to compensate for the sampling.
> 3 votes

---
Tags: rna-seq, statistics, correlation
---


# Question
Title: Difference between paired-end, mate-pair and long read
I writing here because I have some questions for you.

I wondered what the essential differences were between `paired-end`, `mate-pair` and `long read`?

for me `mate-pair` and classic `paired-end` are both `paired-end reads`, with the difference that :

1. For classical `paired-end`: the insert size on classic paired-end is smaller (about 500bp)
2. while the insert size of `mate-pair` is much longer (several Kb) which allows to join the contiguous between them especially is it?

In revenge for the `long-reads`, I imagine that they are simply reads that are synthesized with a large read size but that do not allow like the maite-pair to make junctions between contigs?

# Answer
They are all very different in separate regards, but they all refer to different wet-lab and sequencing protocols/technologies.

First, PE (paired end) reads are typically short (50-300) reads, most often Illumina HiSeq, MiSeq or NovaSeq protocols. Both pairs originate from a single fragment which is sequenced from either end:

```
------> PE1
-------------------------------- Fragment
                         <------ PE2
       ^^^^^^^^^^^^^^^^^^        Insert

```

In contrast mate pairs arise from a fragment that is circularized before sequencing. Now imagine the fragment is circularized:

```
                     -----> MP1
....--------------------------------.......
          MP2 <-----
                    ^Circularized (joined) here

```

Long reads typically refer to PacBio or Oxford nanopore data, with very long reads lengths (up to several hundred kb).
> 4 votes

---
Tags: assembly, read-mapping, long-reads
---


# Question
Title: PacBio long-reads impact in transcriptome de novo assembly?
We are strongly interested in assembly a good transcriptome of reference for a non-model organism and build a local database. We have sequenced the same individual with Illumina (150 millions of pair-end reads) and PacBios IsoSeq v3 (2 SMRT cells, one for shorter transcripts, shorter than 5kb and other for longer transcripts, up to 5kb).

To process long-reads, I have followed the PacBio IsoSeq pipeline proposed in their Github repo (https://github.com/PacificBiosciences/IsoSeq). The final result was removing 70% of the long-reads. Is that normal?

Using this data, I have assembled the transcriptome using only short reads and another combining long- and short-reads. In the end, I have not found any difference... Approx. The same N50, the same number of transcripts assembled, the rate of misassemblies... Does anyone know if PacBio data does not worth for transcriptome de novo assembly?

# Answer
A few comments:

Never use N50 as a metric *especially* for transcriptomes. It has some semblance of relevance for genome assembly, but all that is void for a transcriptome with inherently dynamic lengths.

At the end of your IsoSeq pipeline, you should have (ideally) full-length transcripts. Have you considered that you don't miss much in your IsoSeq data? The biggest gain using IsoSeq is that it should have fewer non-existent isoforms assembled than a De-Bruin graph approach.

Let's say you have the following locus in your genome:

```
5'                                                                 3'
------------=======-----------=======----------=======-------------
Upstream      E1     Intron     E2    Intron     E3      Downstream

```

There are theoretically 7 isoforms:

* E1
* E2
* E3
* E1,E2
* E2,E3
* E1,E3
* E1,E2,E3

A short read assembler has to look at evidence that connects exons to see which isoforms are represented. Let's say you have the following PE reads:

```
5'                                                                 3'
------------=======-----------=======----------=======-------------
               -->..............<--- Read 1
                                 -->...............<-- Read 2
                 --___________-->..<-__________--- Read 3

```

Read 1 and 2 both have a mate in either exon, so that's good evidence that they come from a single transcript. You can also have reads like Read 3, which spans two exon-exon junctions, again good evidence. There are some problems, however, if things get more complicated: if the cumulative length of your exons is greater than the fragment length of your short reads, or you have many exons in a gene, you will never get a *single* read pair that confirms the existence of some transcripts. The same goes for the opposite. If you find a read pair on a single exon, that is not evidence that this single exon is a transcript in and of itself.

Here's where short read assemblers take different approaches. Some make educated guesses about the existence of some transcripts, e.g. by utilizing expression counts in an EM or dynamic programming approach, some just enumerate possibilities.

And that's where IsoSeq comes in. A single, long read from a single transcript helps *a lot* to add evidence and make sure your transcripts actually represent the transcripts that exist in your organism. So, often, your IsoSeq should be *fewer* transcripts than your short read assembly.

In essence, what you should do, is take the IsoSeq and curate those as a set of high-confidence transcripts, i.e. those that you have observed from a single read source. That's not to say that some transcripts that exist only in the short read assembly aren't real, just that you have less definitive evidence for them.

If you are mostly interested in the gene level, that doesn't make a huge difference, but once you need to look into transcript level usage/expression etc it helps a lot to be able to eliminate low confidence transcripts.

Use tools better suited for transcriptome assembly evaluation: Detonate, TransRate and (Busco) and compare your assemblies using those.
> 5 votes

---
Tags: rna-seq, assembly, transcriptome, long-reads
---


# Question
Title: Standard Way to Preprocess Gene Expression?
I am trying to collect gene expression data for the point of fitting gene regulatory networks. My background is primarily in computer science and so I am finding the biological literature a bit difficult to penetrate.

Specifically, I want to fit two regulatory networks for several cancers in TCGA from tumor and normal samples. I have queries using the `R` package `TCGABiolinks` and was able to follow the tutorial here up to the normalization part.

So I have three questions;

* Is there a standard way to normalize gene expression counts? If so, what is it, and is it available as a library in `R`, `python` or some other language?
* If not, can I just center and scale (z-transform) the samples? I am **not** interested in differential expression, only quantifying **expression** from counts.

# Answer
First off - **do not use FPKM for across-sample normalization**. This is such a common issue that I have a prepared text for it:

> Why not to use FPKM

> They have a poor sensitivity when looking at low expressed genes \[1, 2\]. Short transcripts/genes may have falsely over-estimated expression values \[5\]. They may not normalise adequately across samples \[3\]. The reason for this is: FPKM are normalised across libraries using a single variable, namely the “total number of reads”. This library size factor estimation can easily be biased by highly expressed, highly variable genes, unlike the method used by the other approaches, which model the library size factor based on summary statistics extracted from all genes. In other words as stated in \[4\]: “The median ratio method is usually quite robust. The biggest problem for library size correction is not actually the genes with low counts, but the genes with very high counts, which have high variance and disproportionately influence estimators like the total sum. This is why the total sum is a notoriously bad estimator for library size.” Or the same by Lior’s Pachter himself \[6\]. They have a high False Positive Rate in differential expression analysis \[1, 2\]. A lot of the genes/transcripts deemed to be differentially expressed are not, mostly due to the effects of point 1), i.e. high artefactual variability in the expression and point 2) i.e. inadequate library size estimation.

> On the other hand, approaches such as DESeq, DESeq2, vst or voom have adequate between-sample normalisation, have a controlled FDR and can reliably determine differential expression even for lowly expressed genes if there is enough replication or if an adequate log2 fold change cutoff is selected \[7\]. If you really insist on using cufflinks, consider performing upper quartile normalisation, which corrects for the point 1) issue.

> References \[1\] Dillies M.A., Briefings in Bioinf. 2012 \[2\] Soneson & Delorenzi, BMC Bioinf. 2013 \[3\] Lior Pachter’s lecture at CSHL (Cold Spring Harbour Laboratory) from minutes 32 on: https://www.youtube.com/watch?v=5NiFibnbE8o \[4\] From Answer: "Over-correction" in the size-factors of the DESeq2 package: \[5\] http://www.biomedcentral.com/1471-2105/14/370 \[6\] https://liorpachter.wordpress.com/2014/04/30/estimating-number-of-transcripts-from-rna-seq-measurements-and-why-i-believe-in-paywall/ \[7\] http://rnajournal.cshlp.org/content/22/6/839

What I recommend, is that you get the raw count data from TCGA, create a count matrix using the DESeq2 package in R, use the function `varianceStabilizingTransformation()` from DESeq2 to further process your counts and model the networks from there. Your pipeline already mentions `voom`, which is another good alternative, as is `limma`. All of these are available on Bioconductor.
> 2 votes

---
Tags: r
---


# Question
Title: Why is Mann-Whitney U-Test over an ExpressionSet yielding different results than original paper
I am trying to replicate the results of this paper roughly guided by this pipeline. Basically, we are trying to detect a differential expression between mesenchymal and epithelial cells under cisplatin treatment vs control. However, I am facing a problem when trying to calculate the p-value for the Mann-Whitney U-test. This is the part of the paper I am referring to:

> Robust Multichip Average (RMA) normalization was performed at the transcript level on the results from the Affymetrix Human Gene ST1.0 arrays using Affymetrix Power Tool 1.14.4.1 for all 46 sham- or cisplatin-treated ovarian cancer cell lines. The normalized data were subsequently standardized using ComBat71 to remove the batch effect. In this experiment, the cisplatin treatment assay was performed in triplicate on 20 cell lines, while single assays (without replicate) were performed on the remaining 26 cell lines. Taking advantage of the triplicate data, potentially fragile probes with strong variations (an s.d. of \>0.2) within the triplicates were removed, decreasing the probe number from 33 297 to 21 329. To perform a fair comparison, the triplicate data were then log-averaged into one value so that one result for each cell line could be used in the following analyses. To detect differential responses to cisplatin between epithelial- and mesenchymal-like cell lines, the transcriptomic responses to cisplatin were computed by subtracting the gene expression value of control (cisplatin untreated) cells from that of cisplatin-treated cells. Mann–Whitney U-test (P\<0.01 as a cut-off value) was subsequently used to detect the differential transcriptomic responses between the expression changes by cisplatin treatment in epithelial-like cell lines with those in mesenchymal-like cell lines (Supplementary Table 3).

My p-values for a given gene are astronomically different to that of the paper... Here is my code:

expr\_f is the expression values matrix from the Expression set (log-transformed). SDRF\_f is the pData from the same expression set. I am feeding into the wilcox.test function a vector of expression from epithelial cells (epi\_values) and mesencymal cells (mes\_values) for the same probe (8143663).

```
 genelist=rownames(expr_f)
epi_samples = rownames(subset(SDRF_f, SDRF_f[,"Classification"]=="Epithelial-like"))
mes_samples = rownames(subset(SDRF_f, SDRF_f[,"Classification"]=="Mesenchymal-like"))

epi_values = expr_f["8143663", epi_samples]
mes_values = expr_f["8143663", mes_samples]

wilcox.test(mes_values, epi_values, exact=TRUE)

```

In the paper, their p-value for this probe is 0.000372259258165627 while mine is 0.8785. The difference makes no sense at all to me. I don't know if I am not quite getting the concept of the u-test, if I am using the function wrong, or if it's something else. Any input is appreciated

# Answer
> 1 votes
I think you may need to check your pipeline in this regard. I ran the analysis and replicated as closely as I could their M&M (I didn't filter the "fragile probes" since you anyway care only about a specific probe). In any case, I get something relatively close to their result (P=0.004059), which is still one order of magnitude off, but without more detailed M&M, this is probably as close as I can get to their analysis. Here's what I did (if you prefer to read the notebook html, that is here)

library(Biobase)
library(GEOquery)
library(sva) # ComBat
library(oligo)
library(googlesheets4)
library(stringr)
library(dplyr)

#' Get the GEO dataset mostly for the phenoData annotation
gset <- getGEO("GSE47856", GSEMatrix =TRUE, getGPL=FALSE)
adf <- Biobase::phenoData(gset[[1]])
sampleNames(adf)

#' Read RAW cel files from the GEO data supplementary rar archive
#' https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE47856&format=file
#' In my case extracted to ~/Downloads/cel
#' Ensure that we read them in the same order as the phenoData object
celfiles <- sapply(sampleNames(adf), function(f) {
  pat <- paste0(f, ".*.CEL.gz$")
  dir("~/Downloads/cel", pattern = pat, full.names = TRUE)
})
raw_data <- oligo::read.celfiles(celfiles,
                                 phenoData = adf,
                                 experimentData = experimentData(gset[[1]]))

#' Get the raw expression and RMA on transcript clusters
rma_expr <- oligo::rma(raw_data, target = "core")
rma_expr_assay <- assayData(rma_expr)$exprs

#' Use ComBat for batch effect removal
batch_data <- phenoData(rma_expr)[["batch:ch1"]]
cb_expr <- sva::ComBat(rma_expr_assay, batch_data)

#' Get supplementary table 1 from the paper (which I copied on this Google Sheet)
supp_data <- read_sheet("https://docs.google.com/spreadsheets/d/142__1O35zvlttejvz62XQX33uLRJhqIDvb604smVSRo")

#' Match the 'cell line' phenotype from the phenoData to the sample info of
#' the supplementary table and extract the cell classification. Also remove
#' all white space and convert everything to lower case for easier matching.
#' Some cell lines also have a/b attached, so I'm dropping that to match them
#' up. Mayb not be correct to do that here and treat all a/b the same (?)
cell_line <- str_replace(
  str_replace(
  str_remove(
    str_to_lower(adf[["cell line:ch1"]]),
    " +"),
  "^tay$", "taya"),
  "dov13[ab]", "dov13"
)

sample_info <- str_remove(str_to_lower(supp_data$Sample), " +")
classification <- supp_data$Classification[match(cell_line, sample_info)]

#' Get the treament from the phenoData object
treatment <- adf[["treatment:ch1"]]

#' We'll need the geometric mean later (R has no built-in)
gm_mean = function(x, na.rm=TRUE){
  exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x))
}

#' Set up a tibble (==data frame) with all that info for your probe
tibble(Expression = cb_expr["8143663", ],
       Line = cell_line,
       Treatment = treatment,
       Classification = classification) %>%
  # Summarise replicated data with the geometric mean
  group_by(Line, Treatment, Classification) %>%
  summarise(Emean = gm_mean(Expression)) %>%
  # Make sure data is in order
  arrange(Line, Treatment, Classification)%>%
  # Calculate the difference between the cell types
  group_by(Line, Classification) %>%
  summarise(EmeanDiff = diff(Emean)) -> xdata

#' Test the difference
xdata <- split(xdata$EmeanDiff, xdata$Classification)
wilcox.test(x = xdata$`Epithelial-like`, y = xdata$`Mesenchymal-like`,
            exact = TRUE)

#' Session info
sessionInfo()
```

# Answer
> 0 votes
The most likely explanation here is paired versus unpaired.

Mann-Whitney is normally associated with unpaired data, but the authors may have performed a paired test. You performed an unpaired test and thus there is a massive discrepancy in results. The authors should have stated the effect of pairing - it that is what they did - but they are under no strict obligations to do so. You however have performed a Wilcoxan rank sign test, albeit the results shouldn't be so different.

---
Tags: r, microarray
---


# Question
Title: Error in cor(exprs(gse), use = "c") : no complete element pairs
I'm trying to plot a heat map for a "GSE133399", I retrieved the GSE using GEOquery package and then assigned to an object and then tried to plot a heatmap, but I was not successful, I used the following code in RStudio:

```
install.packages("BiocManager")
BiocManager::install("GEOquery")
library(GEOquery)
## change my_id to be the dataset that you want.
my_id <- "GSE133399"
gse <- getGEO(my_id)
length(gse)
gse <- gse[[1]]
gse
pData(gse) ## print the sample information
fData(gse) ## print the gene annotation
exprs(gse) ## print the expression data

## exprs get the expression levels as a data frame and get the distribution
summary(exprs(gse))
exprs(gse) <- log2(exprs(gse))
boxplot(exprs(gse),outline=FALSE)
library(dplyr)
sampleInfo <- pData(gse)
sampleInfo
## source_name_ch1 and characteristics_ch1.1 seem to contain factors we might need for the analysis. Let's pick just those columns

sampleInfo <- select(sampleInfo, source_name_ch1,characteristics_ch1.1)

## Optionally, rename to more convenient column names
sampleInfo <- rename(sampleInfo,group = source_name_ch1, patient=characteristics_ch1.1)
sampleInfo
install.packages("pheatmap")
library(pheatmap)
## argument use="c" stops an error if there are any missing data points

corMatrix <- cor(exprs(gse),use="c")

```

I get this error:

```
Error in cor(exprs(gse), use = "c") : no complete element pairs

```

I also tried this:

```
corMatrix <- cor(exprs(gse),use="c", use="complete.obs")

```

I get this error:

```
Error in cor(exprs(gse), use = "c", use = "complete.obs") :
  formal argument "use" matched by multiple actual arguments

```

I want to achieve this plot:

# Answer
> 2 votes
The GSE133399 dataset contains both ATAC-seq and RNA-seq data, none of which can be queried with GEOquery which works only for microarray data.

`exprs(gse) ## print the expression data` returns an empty object as you might have noticed, therefore the errors in the lines that follow this command.

On the bottom of https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE133399 at the `Supplementary file` section you can download FPKM-transformed count data and use these. I would be careful though as FPKM is not ideal for inter-sample comparison. I personally always (if feasable) download the raw data and then process myself, but for a guess on what is going on the FPKMs probably will do.

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


# Question
Title: Finding motifs: fasta file with 10,000 sequences
I am new to Python. I am trying to parse a fasta file containing 10,000 sequences to look for motifs (microsatellites in particular). I tried using Seq Utils to parse my sequences for a particular motif, that is, (TG)20. I suspect because I used the string function (str), it treated the entire fasta file as a one single fasta sequence. I then tried using BruteForce on the motif (TG)17 and found that an error appeared which said that "...the object of type 'FastaIterator' has no len()...". The script however worked when only a single sequence was in the fasta file.

In essence, I need some help figuring out how to:

* Prepare an instances object(is this correct?) with microsatellite motifs (for starters, I intend to parse the sequences for dinucleotide repeats (TG)10,(AC)10,(AT)10 and (AG)10;
* Parse my 10,000 long sequence file using this object(?); and
* Print an output specifying which of my sequences contain microsatellite motifs.

Appreciate any help I can get on this.

#Seq Utils to find Motifs

```
>>> pattern = Seq("TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG")
>>> sequence = Seq(a)
>>> results = SeqUtils.nt_search(str(a),pattern)
>>> print(results)
['TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG', 38, 40, 42, 44, 46, 48, 50, 52, 54,
 56, 58, 68682, 68684, 68686, 370986, 370988, 370990, 370992, 370994, 370996, 43
8620, 438622, 438624, 438626, 438628, 438630, 438632, 438634, 438636, 438638, 43
8640, 556703, 556705, 556707, 784973, 784.... # it's too long so I am just pasting the first few lines

```

#BrutForce to find Motifs in a fasta file containing 10,000 sequences

```
>>> def BruteForce(s, t):
...     occurrences = []
...     for i in range(len(s)-len(t)+1):
...         match = True
...         for j in range(len(t)):
...             if s[i+j] != t[j]:
...                 match = False
...                 break
...         if match:
...             occurrences.append(i)
...
...     print(occurrences)
...
...

>>> t = "TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG"
>>> len(t)
34
>>> s = SeqIO.parse("group23.fas", "fasta")

>>> BruteForce(s,t)
Traceback (most recent call last):
  File "<console>", line 1, in <module>
  File "<console>", line 3, in BruteForce
TypeError: object of type 'FastaIterator' has no len()

```

#BrutForce to find Motifs in a fasta file containing a single sequence

```
>>> s = SeqIO.read("group23v.fas", "fasta")
>>> BruteForce(s,t)
[0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 54, 56, 58, 60, 62, 64, 66, 68, 70, 79, 98, 117, 175, 193, 198, 217, 221, 277, 279, 282, 289, 294, 304, 311, 335, 381, 437, 464, 468, 510, 514, 522, 593, 602, 622, 634, 647, 658, 678, 704, 727, 750, 773, 804, 889, 895, 908, 911, 922, 937, 969, 975, 978, 1002]

```

# Answer
`Bio.SeqIO.parse()` returns a `SeqRecord iterator`. To access the sequence of this sequence record, you need to use the `.seq` attribute so you should update your `s` with `s.seq`.
> 1 votes

# Answer
To expand on haci's answer (and since you mentioned being new to Python) you can loop over the iterator object and use the `seq` attribute like this:

```
pattern = Seq("TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG")
for record in SeqIO.parse("group23.fas", "fasta"):
    results = SeqUtils.nt_search(str(record.seq), pattern)
    print(results)

```

(More info on iterators here.) That will work the same whatever the number of sequences in group23.fas.

I'd be cautious about looking for microsatellites with fixed strings. If the sequences you're looking for are short and have imperfections among the repeats, you might miss them, depending on how many perfect repeats your search pattern has. For example with these stub example sequences:

```
>seq1
TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG
>seq2
ACGATGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG
>seq3
TGTGTGTGTGTGTGTGTGTGTCTGTGTGTGTGTGTGTGTG

```

The nt\_search output gives you:

```
['TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG', 0]
['TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG', 4]
['TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG']

```

That is, one match at position 0 of the first sequence, one at position 4 of the second, and no matches for the third (there's a C I put in the middle). This might not be relevant for you if you're looking microsatellites much longer than your search patterns, but for shorter ones it might be a problem.
> 1 votes

# Answer
Here's a solution that I came up with after help from both haci and Jesse. It involves using the following modules: os, glob, SeqIO, Seq, SeqUtils.

```
import os, glob

os.chdir("C:\python38\Lib\site-packages\Bio\SeqIO")

from Bio import SeqIO

from Bio.Seq import Seq

from Bio import SeqUtils

patttern = Seq("ACACACACACACACACACAC")

folder_path = 'C:\python38\Lib\site-packages\Bio\SeqIO' #specified the folder where my sequences were stored

fasta_paths = glob.glob(os.path.join(folder_path, '*.fasta')) #glob module finds all the pathnames matching a specified pattern i.e. '*.fasta'

with open("motif4.txt", "w") as f:
... for fasta_path in fasta_paths:
...     print(fasta_paths)
...     for seq_record in SeqIO.parse(fasta_path, "fasta"):
...         results = SeqUtils.nt_search(str(seq_record.seq),pattern)
...         f.write(str(results)) #replace with seq_record.seq to save the list of node names of all sequences that were parsed
...         f.write("\n") #ensures results of the parse through all the sequences in the fasta_path is written to a new line
... f.close()

```

The nt\_search output gave the following output (similar to what Jesse showed):

```
1 ['ACACACACACACACACACAC']
2 ['ACACACACACACACACACAC']
3 ['ACACACACACACACACACAC']
4 ['ACACACACACACACACACAC']
5 ['ACACACACACACACACACAC']
6 ['ACACACACACACACACACAC', 16636, 16638, 16640, 16642, 16644, 16646, 16648, 16650, 16652, 16654, 16656, 16658, 16660, 16662, 16664]
7 ['ACACACACACACACACACAC']
8 ['ACACACACACACACACACAC']
9 ['ACACACACACACACACACAC']
10 ['ACACACACACACACACACAC']
11 ['ACACACACACACACACACAC', 9174]

```

It would miss sequences that did not match (AC)10 or those with imperfect repeat motifs. However, in cases where the motif was \> (AC)10, nt\_search listed down all locations where the motif was found sequentially. For example, in line 6, all locations were listed and showed the motif was in fact 36 nucleotides long in total (pos #16,636 to #16,673).

You can then import the txt. files into a databasing app and analyse the results.
> 1 votes

---
Tags: python, sequence-analysis, motifs
---


# Question
Title: Add tags to the read pairs in a bam file
I have a bam file from paired-end sequencing, in which each R1 has the `RX` tag assigned (via `UMItools group`) denoting a Unique Molecular Identifier (UMI).

I want to call molecular consensus reads using `fgbio` but `fgbio GroupReadsByUmi` gives me an error, because it expects both R1 and R2 to have the `RX` tag set.

I would like a tool that will go through my bam file and assign each R2 the same `RX` tag as the R1.

# Answer
> 1 votes
The best thing would be to figure out how to do this in UMItools. I don't know that tool at all, but it seems that it is supposed to be able to deal with / tag both reads in paired data. e.g. did you use `--paired` flag, etc.

I would suggest investigating that route first before you look at the next part of my answer.

---

Per Devon Ryan's comment suggestion I would suggest a workflow similar to this:

* sort bam file by query name (`samtools sort -n`) to get the R1s immediately following the R2s
* Write a short script/command that
    1. reads through the bam file and for each line
    2. looks up the tag in the R1, prints R1 line
    3. reads the R2, substitutes the RX tag in, prints R2 line
* pipe that through `samtools view` and possibly `samtools sort` back to a bam file of the same type that you started with.

My pysam is a little rusty but something like this might work

```
import pysam

bam = pysam.AlignmentFile("bam", "rb")
# print header
print(bam.header)

tag = None
for read in bam:
    r2 = False
    if tag:
        read.assign_umi_tag() # assign to read- not sure where/how
        r2 = True
    else:
        tag = read.lookup_umi_tag()  # i don't know where this tag is in bam format
    print(read)
    if r2:
        tag = None

```

I haven't tested that in any way, caveat emptor.

I am not personally familiar with UMItools and how it assigns tags so can't suggest anything more specific, sorry.

---
Tags: sam, umi
---