title: Reading in external single cell data

text:
# Question
Title: Reading in external single cell data
I'm trying to read in an external single cell dataset from https://www.nature.com/articles/s41467-020-16164-1, but I am having trouble reading in the count matrix.

counts found here: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE154778

This is how I've been reading in the raw counts.

```
cts = read.table("GSE154778_RAW.tar", sep=',', header=T, stringsAsFactors=F, check.names = FALSE)

```

The output is this:

```
> head(test)
                                                                                GSM4679532_K16733_barcodes.tsv.gz
1                                                                                xVZ\xbe\xc6pܧj\x95I\xed\030\xee\xfe\xf8ޔ\x93\xa7\xbc\xf2\023\xed{\x8e\x98)w\\\xef:ܣN\xdfrr\xc7\xf5\f\xaa\xcb\001\xfa[\002\021@|a\x91\024\xab\xc5^\017\xc0O\xbd`7\x97&m\xb9\x81\xe8R\x97{\x9f\023@\x9a\xf2\x81I\0368\x80N\027\xfa\xf8\x89fǾ6V\x98h\x9aF\023\x80\xef\x8f\xf2;f\xb1%@\xffif\004(4\xf4\037\020\xaf\x8d\032\037Y\021 1\xdd\xd9V\xc0\xde{\xf0U\017<?\xd6\xe0D\x9b\x96\x8e\xc1\xd4s\xdb\024
2                                                                                *\x9b\xee5\xbb\xb3g\xab%\xa5r\xd7?\x82\x89V\xe5\177\xd5L4g\xe7?\xa2\x96\x8b\xaa\xf0\0215D\xd4\xe1\x87~D\xd9ou\xda\xe7?~\xab
3                                                                                e\xb1\xe9ky\xb0w0Qy\xd7\xdd\xc6uaZ\xcb8\025-N\x96\xb6Y\xac\0015h\xeb\xe5\nL\xeb\xc4\177\x9cX\xbdU\016\xde\xcc\xdeAE\xa7\\w\xfcE_\001\xbb\026\xba\xe7\xf7ܳS'\x96\u038bI\xfe\x9a\xef\xfdm\xa2\xbf_瘺a\030\037\xec\xa3\032\x97\xa2\xea\x9e\xf3\xf40\x8d\xfb#\n\xa9\xb8vE\x94\x9b:\xef\xcb\xe5\xd5\xc69n\xe1\xbd\xcab\025\xc7\xfb|\fT\xb5\xfa\xea4lK\xb2\037)\xdfL\xfc\x8f\xf5yz\xa6D}\001/\xaak?̑\xe9\003\xf1\031\x94\xed~\xfdӳ\xc7B]\xca/\\\xfb\xa3%n
4 \xef\xaa\xc4B\xe0]\xe5\xac6\x9e\xea\xbbe\xf9\x8a\xecch/nV5\xd1p_O\xb6N\034\xa9\xf3x\034\x98&\025\xd75\xa7b\021\xb3v\xd6P\xf7D4~u\xb3\xab(\xf9\x9d\xb29^,\xd8\xf3v`B)\xc3c|>'\xc36Z\016\xb4\0257x8\u07b5fb\xd1(\xb9\xbe\xe7WT.x\xaa5\xe4\xc1\xd6\xe6ny\f֟\fYST.\\\xd0\xf2h\xf6\x8d\xc5Q\x90o\xbc\xadG\004\xb1\x9el\xcb\xd8\037\xef\xd4\024Y\xf1Q.\xe3\x81FGZ\016\026t\xa4\xcaY\xf3\xb7\x95\017,\au\xe8\037\xcf\a\xb3\ab\001\x8d.\aճ\x8f\xe6\xc4h\xe5n8\xf7Sdg\xc5Z\xef\xaa\xcb!\xd6\xd0ӳ\x92\xa2\xc4Wݥ?\xcc\xc4\016\031\023\xaa\xe7\xb6-H^X\005\xb1ӡ\x9b\xe6\xfb\xf0\x9b\x9d\027q\xcbq\x9aEm%\036\xb3\xe9$f\xc8\xd6\005\xdb\xee\x8b\xec\xf5\xb9\xef'\xaa~\xad\xb3\xf9\x8fg\xb3\025\xe5h\xc9\xe9_\xee7v\x87\xa3\xf5\177,ǌ\x9c\016Z\f\xe6\xf353\027\0367=\xc0\002{\x89\037,\017\xd3Œ\030\034D/\xeeb\x8ct\033\xe5UD\xedϧ\xb9\xf1\xbe\xe2\034L\xa7o\xea\xff\xb9\xecŭ\x9b\036Õ\xdbM+>Su$\\\xe0\xda2\xdbl?\xa2\xcc\v\xf4rS\xe4\027\x8aEdR9f\xb9N\022\xb4И@,e\xfe<\x95gN\xb4Ć\x9b\xc7\xeb\032\xc7P~\xa2\034\xd8]\x98\032\te\xff\xd8!vr\xb0\xa6\x92K\xf0\xfd`*\xb3\xe8\aq\x90\x9e'\xe6FI\x84\xa3<B\xeddE\xefpg\001V[\xbc\xb7u\xb04\003\xe3K\xb6h\xc0ild`\x92\025W\xe5#J\177\xa6\x95\xabZ.\x96\xec\016w\x93d÷+]Di\xbaJ\xa2\003\xf5\021m\xa4\x80J\022\f\027QqE\021r\x94\xf0\xbd\037\xd6ye\x970Fm\xa1ˎO\035\xfe\xe5\xb8\xf2/\x89\x9bX$\xa9b\026\037B~\xe2\023k\033\xbaB\xbcɪ\xa7\xe4f\xf2\027\x8f\xa2t\xa1\xee\xbd\xec\x80\xe5;\037\xbc\xfeC\x94\xed\x91UQ\x81~\xbe^N<\003&\xceDJE\xab\031\xa7\xf9\x9a\xa2\xf2%\024}\x9e.g\x9a0J&lj\xd2\032\xc55~\022\xa3ޭ\xd3|\xd1,k#\017ֶ\xfbg\026CL\016\x9bx\023Y\x89\xffI\x85\x8c\x89>U=\024|í\xb4\xf2\006\xfa#ʎ\032v\xc6Hps\034\fI\xef`E\xb1\xec\xcb\xf3\xc6\xf2E\xe5ӹ\xee1σc\x8b\x89\xe5\xf2\xd0\xe3\xf2\xeb\xac\xe5\xfb,&N\002\x9b\xe14i\xc5\024y\xc8^Ď\xe73\n\006\022\x94TM\xc5p\xbf(\xb8w~ES\x8d&N\xafr\xefL\f\xbf\x816~$\x99\021A\xfe7m\tƘhX\xe0\xa2\xc2>&\xe5s\xa0\x99\t\004\x88\u07bd\xa34\xf6q\xe4\xf1ٕ\x9b\xea\003\xd7frq(~\xe7\024\035\xcf
5                                                                                \x91\xb7U
6                                                                                \x9f-\xcf\025\xa7\xeb6<\x89\xf55\xc6+\xb9K\xb5\xbaa-[.\xda \xca]<M\xea\xbd%

```

Any suggestions?

Update:

I unzipped the file in HPC but now I'm having trouble reading in the counts as it's all the counts per sample.

```
Data_Processing.nb.html            GSM4679535_T3_genes.tsv.gz     GSM4679540_T9_barcodes.tsv.gz      GSM4679544_Y00014_matrix.mtx.gz
Data_Processing.Rmd                GSM4679535_T3_matrix.mtx.gz    GSM4679540_T9_genes.tsv.gz         GSM4679545_Y00016_barcodes.tsv.gz
GSE154778_dgeMtx.csv.gz            GSM4679536_T4_barcodes.tsv.gz  GSM4679540_T9_matrix.mtx.gz        GSM4679545_Y00016_genes.tsv.gz
GSE154778_RAW.tar                  GSM4679536_T4_genes.tsv.gz     GSM4679541_T10_barcodes.tsv.gz     GSM4679545_Y00016_matrix.mtx.gz
GSM4679532_K16733_barcodes.tsv.gz  GSM4679536_T4_matrix.mtx.gz    GSM4679541_T10_genes.tsv.gz        GSM4679546_Y00019_barcodes.tsv.gz
GSM4679532_K16733_features.tsv.gz  GSM4679537_T5_barcodes.tsv.gz  GSM4679541_T10_matrix.mtx.gz       GSM4679546_Y00019_genes.tsv.gz
GSM4679532_K16733_matrix.mtx.gz    GSM4679537_T5_genes.tsv.gz     GSM4679542_Y00008_barcodes.tsv.gz  GSM4679546_Y00019_matrix.mtx.gz
GSM4679533_Y00006_barcodes.tsv.gz  GSM4679537_T5_matrix.mtx.gz    GSM4679542_Y00008_genes.tsv.gz     GSM4679547_Y00027_barcodes.tsv.gz
GSM4679533_Y00006_genes.tsv.gz     GSM4679538_T6_barcodes.tsv.gz  GSM4679542_Y00008_matrix.mtx.gz    GSM4679547_Y00027_genes.tsv.gz
GSM4679533_Y00006_matrix.mtx.gz    GSM4679538_T6_genes.tsv.gz     GSM4679543_Y00013_barcodes.tsv.gz  GSM4679547_Y00027_matrix.mtx.gz
GSM4679534_T2_barcodes.tsv.gz      GSM4679538_T6_matrix.mtx.gz    GSM4679543_Y00013_genes.tsv.gz     Lin_et_al_GMedicine_metadata.csv
GSM4679534_T2_genes.tsv.gz         GSM4679539_T8_barcodes.tsv.gz  GSM4679543_Y00013_matrix.mtx.gz
GSM4679534_T2_matrix.mtx.gz        GSM4679539_T8_genes.tsv.gz     GSM4679544_Y00014_barcodes.tsv.gz
GSM4679535_T3_barcodes.tsv.gz      GSM4679539_T8_matrix.mtx.gz    GSM4679544_Y00014_genes.tsv.gz

```

I've created this loop to read through all the samples.

```
counts_list <- list()
for (sample in samples){
    print(sample)
    filename_header <- sample
    barcode_file <- file.path(wd3, paste(filename_header, "_barcodes.tsv", sep=""))
    gene_file <- file.path(wd3, paste(filename_header, "_genes.tsv", sep=""))
    matrix_file <- file.path(wd3, paste(filename_header, "_matrix.mtx", sep=""))
    counts <- Matrix::readMM(matrix_file)
    barcodes <- read.csv(barcode_file, sep="\t", header=F)
    barcodes <- barcodes$V1
    colnames(counts) <- barcodes
    genes <- read.csv(gene_file, sep="\t", header=F)
    genes <- genes[,-3]
    colnames(genes) <- c("ENSG", "Name")
    rownames(counts) <- genes$ENSG
    genes <- genes[!duplicated(genes$Name),]
    counts <- counts[genes$ENSG,]
    rownames(counts) <- genes$Name
    counts_list[[sample]] <- counts
}

```

Receive this error when I run:

```
counts <- do.call(cbind, counts_list)
Error: Matrices must have same number of rows in cbind2(argl[[i]], r)

```

I added:

```
genes <- genes[,-3]

```

to ensure that all the gene file has only 2 columns as I noticed that one sample had 3 variables but I'm still erroring.

# Answer
Tar is an archive format so as you figured out by now you have to unpack it.

The files you now see are as follows:

* `mtx.gz` are compressed count matrices which you load into R using the `Matrix::readMM()` function
* `genes.tsv.gz` are the associated rowwise data so the "genes" which you can load with `read.delim` and then use `rownames()` to assign them to the matrix generated by `readMM()`
* `barcodes.tsv.gz` are the associated column data so the samples or "cells", use `colnames()` to assign them to the matrix

I would use `list.files` to create a vector of files and then either a loop or `lapply` to iterate over them for the loading. Does that make sense to you?

**Edit after OP received this dimensions error**: This then means that these files come from different experiments and have different upstream processing. THere is nothing you can do about it but figuring out which files come from which toplevel experiment and then only combine those of one experiment, respectively. The procedure towards reading data is as I lined out above.
> 3 votes

# Answer
It looks like it is 10X data and I strongly suggest a dedicated single cell analysis package such as Seurat. You will most likely end up using one of these anyway for data analysis.

Read10X() can be a good start. I don't remember whether it requires dedicated folders per sample though. Even if this is the case, you can create individual sample folders with a simple bash script, can be done within R as well.
> 1 votes

---
Tags: r, seurat, single-cell, data-preprocessing
---


# Question
Title: What are some ways to check metagenomic bin quality?
I am new to metagenomic binning. I've used CheckM in order to estimate completeness and contamination values, and most of my bins of interest appear to have good values. My workflow was pretty straightforward (BBDuk for quality control \> Megahit for assembly \> vamb for binning). What else could I do / where do I go from here? (especially regarding quality control)

# Answer
Some papers check 5S, 16S, 23S and ≥18 tRNA genes for completeness. For example, you can see this paper.
> 1 votes

# Answer
There are many other tools. If you want to get nice statistics and graphs I recommend Amber, there is also Anvio for manual refinement of the bins.
> 1 votes

---
Tags: assembly, metagenome, quality-control
---


title: Offline Interactive Genome Map

text:
# Question
Title: Offline Interactive Genome Map
I use the CGView Server for constructing circular genome maps. It is really a convenient tool for me to explore circular genome maps because it's interactive. I can search a certain region and it can automatically hover to that region. However, the job running on CGView Server will be deleted after 7 days. Are there any alternatives to CGView Server? Preferably something that's offline so I can save it for longer than 7 days. I've checked the website on CGView, but I think the java package for CGView is not interactive. Thanks in advance!

# Answer
SnapGene Viewer: https://www.snapgene.com/features/

MacVector: https://macvector.com/

ApE: https://jorgensen.biology.utah.edu/wayned/ape/
> 1 votes

# Answer
I've been using Unipro Ugene and Snapgene Viewer (not for circular genomes, but they both have the option).

Snapgene Viewer is more user-friendly, but I think Ugene might have more tools that are free (if you want to copy-paste pieces of DNA, for instance, you can't do that in Snapgene Viewer, it would require the paid or demo version of Snapgene)
> 1 votes

---
Tags: genome
---


title: Nanopore variant calling

text:
# Question
Title: Nanopore variant calling
So far I haven't done any variant calling as such. Nanopore I have used for 16s microbiome data.

Now My question/doubt so how do I proceed for nano-pore virus sequencing data

Steps:

* I get fast5 files
* The i perform base-calling using guppy
* Get fastq files

Now after doing the above is it possible to use variant calling that is used for Ilumina once I get the fastq files?

# Answer
Tools for illumina are not really suitable for long read sequencing. For nanopore data I would suggest LongShot, Clair3 or DeepVariant-PEPPER.
> 1 votes

# Answer
Sure it is possible, but I guess it won't be more reliable than a method specifically designed for Nanopore. But if you use Nanopore I guess accuracy is not your main concern.
> -1 votes

---
Tags: variant-calling, nanopore
---


title: BLAST server download

text:
# Question
Title: BLAST server download
I discovered that Phytozome's BLAST results look quite nice (please see below). By any chance, does anyone know where to download this BLAST server in order to install it locally for organisms which are not part of Phytozome?

Thank you in advance,

# Answer
I don't know about the BLAST server you're referring to unfortunately.

But another BLAST server, SequenceServer looks quite nice too. And has a bunch of helpful visualisations that you can download directly in png or svg. Source code and installation instructions for this one are at http://sequenceserver.com

My opinion is biased because my lab makes SequenceServer. But there is also social proof: dozens and dozens of community websites use it, and it get cited regularly.

Apologies for spamming a bunch of screenshots. I think they're helpful to understand what SequenceServer can do & I'm proud of what we've built.

(it's open source so you can install it on your own server, or use our cloud BLAST hosting service).

## Overview of how hits align to query

## Lengths of queries and hits

## Kablammo pairwise comparison of hit-query alignment

## Circos overview of BLAST results

## Search interface (automatically detects query type & chooses BLAST algorithm)

## Example BLAST result overview
> 1 votes

---
Tags: sequence-alignment, blast, webservice
---


title: VCF merge containing CNV

text:
# Question
Title: VCF merge containing CNV
How do I merge VCFs files containing CNVs?

I use vcf-merge, a VCFtools function, and after bgzip and tabix, SAMtools, to index and tab separate variants, but I don't know if it is the right way.

# Answer
> 2 votes
`vcf-merge` (or better: `bcftools merge`) does a column-based merge of VCF files, basically taking on an additional genotype column to each variant. While this sort-of works, it doesn't take into account all the information about variation. There might be low coverage at a region in one VCF file that indicates a deletion, but the deletion is only observable when considering other samples with high coverage at that same region. It's also possible that INDELs could be defined differently in different samples, requiring an additional `bcftools norm` step for left-normalisation prior to the merge.

Assuming this column-based merge is what you want (and if not, then `bcftools concat` is your friend for row-based merging), then it's better to regenerate the merged VCF file as a multi-column file directly. One way of doing this is feeding *all* BAM files into `samtools mpileup` to generate genotype likelihoods, then feeding the result of that into `bcftools call`. Something like the following:

```
samtools mpileup -uf ref.fa in1.bam in2.bam in3.bam | \
  bcftools call -mv -O z -o variants_ref_1-3.vcf.gz

```

# Answer
> 0 votes
depending on your goals you could possibly use tools like https://github.com/mkirsche/Jasmine or simple bedtools intersect on the vcf entries could help (e.g. harmonize after the calling)

or if you have a cohort maybe possibly doing joint-CNV calling or "genotyping CNV" in a cohort is useful http://penncnv.openbioinformatics.org/en/latest/user-guide/joint/ or https://github.com/abyzovlab/CNVpytor#merging-calls-form-multiple-samples

---
Tags: vcf, cnv
---


# Question
Title: Calculating dn/ds ratio between two model organisms
I am trying to calculate the dn/ds ratio between the human genes and c.elegans. However, I am getting NA values in the output file

I used the following script

```
library(biomaRt)

## This gives the full list of all marts, including archives
archives <- listEnsemblArchives()

## For the c.elegans data
ensemblcelegans=useMart(host = 'jan2020.archive.ensembl.org', biomart =
'ENSEMBL_MART_ENSEMBL',
                 dataset = "celegans_gene_ensembl")

filters = listFilters(ensemblcelegans)
attributes = listAttributes(ensemblcelegans)

Celegans_GFList <- scan("c.elegans.tsv", what = character())
Celegans_GFListValues<-getBM(attributes = c('ensembl_gene_id',
'hsapiens_homolog_ensembl_gene',  'hsapiens_homolog_dn',
                                    'hsapiens_homolog_ds', 'hsapiens_homolog_orthology_type'),
                     filters = 'ensembl_gene_id',
                     values = Celegans_GFList,
                     mart = ensemblcelegans)

Celegans_GFListValues1 <- subset(Celegans_GFListValues,
Celegans_GFListValues$hsapiens_homolog_orthology_type == 'ortholog_one2one')
Celegans_GFListValues1$dnds_ratios <-
Celegans_GFListValues1$hsapiens_homolog_dn/Celegans_GFListValues1$hsapiens_homolog_ds
write.csv(Celegans_GFListValues1, file = "Celegans_GFListValues.csv")

```

# Answer
> 0 votes
Ensembl discontinued our calculation of dN/dS scores for orthologue pairs in release 100 (https://apr2020.archive.ensembl.org/index.html). However, even when we did perform these calculations in Ensembl 99 and earlier, we only did this for pairs of orthologues that diverged within the last ∼100 million years, namely between pairs of mammalian, avian or tetraodontiforme genomes. More information can be found here:

https://academic.oup.com/database/article/doi/10.1093/database/bav096/2630091

---
Tags: r, biomart
---


title: Understand this Kmer plot from Merqury?

text:
# Question
Title: Understand this Kmer plot from Merqury?
The attached figure is generated based on Illumina reads from multiple individuals compared to genome assembly. Looks like there are a lot of kmers are reads only (grey colored). Also, a blue peak (2x kmers), does it represents high heterozygosity? GenomeScope gives 1.5% of heterozygosity result. But it seems the blue peak normally should be near the 70x on the x-axis. Thanks.

# Answer
> 1 votes
Assuming a diploid genome, usually, we aim to build a monoploid reference, representing a single haplotype.

Practically speaking, we need to include only one of the two k-mers that are representing heterozygous sites, however, we would desire all the homozygous k-mers to be included exactly once.

Looking at your plot, I would guess 1n coverage is somewhere around 18x, and 2n about 36x. If that's right, your plot is suggesting that your assembly contains lots of uncollapsed haplotypes, which is quite understandable given very high heterozygosity you estimated with genomescope. This is also supported by many homozygous k-mers are there twice (the blue bump). These are probably regions in between heterozygous sites, or flanking them.

Anyway, I would recommend something like purge haplotigs to collapse haplotypes into single monoploid reference and redo the plot. The blue bump should disappear and the black and red distributions should have about the same height around 18x, while nearly all the 36x k-mers should be red (once).

Also, if the interpretation of the plot still would not be clear, you can check this part of a lecture by Lucía Campos-Dominguez on the topic (it's a different software called KAT she talks about, but it is using the same principle as Mercury).

---
Tags: genome, assembly, illumina, k-mer
---


# Question
Title: Pymol python script for selecting a residue's sidechain and calculating its center of mass
I trying to write a pythonic PyMol script that loads a PDB, cmd.iterates through the PDB and stores the residues (three letter code) into an array. I then want to loop through this array and calculate the center of masses (COM) of each aa's sidechains (except gly) and add it into a second array. I am able to calculate the COM of a residue as a whole, i.e. including backbone atoms. But what I want to calculate is the COM of the sidechain. Please point me in the right direction.

```
path = 'Path to PDB dir'
for filename in glob.glob(os.path.join(path, '*.PDB')):
    cmd.load(filename)
    stored_residues = []
    cmd.iterate('name ca', 'stored_residues.append(resn)')

    stored_residues_COM = []

    for i in stored_residues:
         stored_residues_sidechain_COM.append(centerOfMass('resn %s' % i))

```

# Answer
I am guessing you used the command `pymol.cmd.centerofmass`? If so the selector is the ticket to doing it by selecting everything of the residue bar the backbone atoms, e.g. f'resi {r} not name N+CA+C+O+H'

Parenthetically, the `cmd.iterate` command is really ugly and counter-intuitive —it's heavily used because it is older. The command `get_model(selection)` in pymol2 module is so much nicer.

```
for atom in pymol.cmd.get_model('name CA').atom:
     atom.resi #note it is type str due to insertion codes!
     atom.resn
     atom.coords # List[float]

```

Second parenthesis. A `and not element H` is not needed or beneficial to `pymol.cmd.centerofmass` as this takes account of mass and occupancy. There is no need to select only heavy atoms as a proton weighs nothing. The barycentre and centroid are the same in geometry where there is no mass to the points, but here mass is present.
> 1 votes

---
Tags: python, pymol
---


title: STAR aligner multiple fastq files

text:
# Question
Title: STAR aligner multiple fastq files
I’m using STAR to align fastq files from SMART-seq2. I have raw data folder containing sub-folders with samples names the sub-folders each contain fastq file. How can I make a bash command in order for STAR to go in each of sub-folders and make alignment for the fastq file and retrieve the bam files in one folder?

Folder structure:

```
/plates/plate_1/rawdata/A1/A1_R1.fastq
/plates/plate_1/rawdata/B1/B1_R1.fastq
...
/plates/plate_1/rawdata/O24/O24_R1.fastq
/plates/plate_1/rawdata/P24/P24_R1.fastq

```

It’s 384 sample folders A1, A2, B3, F5,,, so on

# Answer
It seems the sequencing is single-end, and there is only one fastq per sample. If so:

```
staridx=/path/to/star/genome/index

mkdir /plates/star/
cd /plates/star/

files=$(find /plates/plate_1/rawdata/ -name "*_R1.fastq")

for i in ${files}
do

  sample=$(basename $i "_R1.fastq")

  STAR --genomeDir $staridx --runThreadN 8 --outFileNamePrefix ${sample}. \
    --readFilesIn $i --outSAMtype BAM

done

```

All fastq files will be mapped individually, and all bam files will be saved to `/plates/star/`.
> 2 votes

---
Tags: sequence-alignment, genome, bash, star, loop
---


title: Nanopore FASTQ header specifications

text:
# Question
Title: Nanopore FASTQ header specifications
Is there a specification anywhere for how a Nanopore FASTQ header should be formatted?

I am creating strategies for generating test FASTQ files and would like to ensure I am generating "correct" nanopore sequence identifiers.

# Answer
> 3 votes
The headers are defined by the basecaller (guppy, Albacore, bonito, etc) during basecalling. They have changed many times over the years, but generally they start with a readID of `@+UUIDv4` followed by multiple entries separated by spaces, of form key=value

For example:

```
@945da76d-b7ef-4c65-91a1-e69c85832185 runid=2cf9e7c4dd1888d9a1090ebc394dd5b7cb2fe3f0 read=126 ch=181 start_time=2017-10-17T06:43:19Z barcode=barcode01

```

This is still valid, however, there is some movement to add sam tags to parse data through to sam files when aligned.

---
Tags: nanopore, fastq
---


# Question
Title: How to extract a subset of individuals from a .bgen file?
I have 10,000 controls and 3000 cases and the .bgen file has genotype of 200,000 individuals. I would like to extract only those individuals from the .bgen file that are either cases or controls using plink for downstream analyses. How can i do that?

# Answer
> 2 votes
You can use qctool for this. If you create a list of IDs in `list.txt` (one ID per line), you can filter the .bgen file(s) by running

```
 qctool -g yourdata.bgen -s yourdata.sample -incl-samples list.txt -og yourdata_filtered.bgen

```

---
Tags: file-formats, plink, filtering
---


# Question
Title: Chou-Fasman Method Overlapping Regions Problem - Secondary Structure Prediction
According to the Chou-Fasman method, when there is an overlapping region on the amino acid sequence, we compare the scores of alpha-helix and beta-sheet and chose the highest one. Let's say I have this amino acid sequence: **ABC**.

As a result of computations, **AB** is identified as alpha-helix, and **BC** is identified as beta-sheet. This means, I have an overlapping region, which is **B** in this case. The algorithm says that we need to compare the scores of the overlapping region, and select the secondary structure with the higher score, let's say alpha-helix in this case.

Would the resulting secondary structure be $\alpha$$\alpha$$\beta$ or $\alpha$$\alpha$*_*?

In other words, do we cancel out the secondary structure identifications when we encounter overlapping regions?

Thanks in advance.

# Answer
> 1 votes
Sure, an isolated residue has a ψ and φ angles, so you could technically have an isolated residue in a given SS conformation, but that does not mean much if its neighbours' backbones aren't in that same conformation. Also, sure, in Rosetta design algorithms residues are thrown into ABEGO bins based on ψ and φ angles alone and as a result one can get isolated α-helix bin residues (in what is a loop): but that's a design algorithm that will want to correct this.

In your particular example, IIRC that algorithm uses a window of 4 or 6 residues for SS prediction. So it will surely be (αα)αα\_. But I am pretty certain you cannot have a ααααββββ switch without a loop or turn in nature —I might even hazard a guess that for design it is not Koga-compliant even if not specifically addressed. And with the Chou-Fasman algorithm you'd need to have the windows spanning the switch score above the threshold in the correct way. I am sure there's a combination that allows it as after all we're talking of an ancient 70s prediction method based on probabilities alone...

---
Tags: protein-structure, secondary-structure
---


# Question
Title: Searching for adapter sequences in FASTQ files - metgenomics
I have recently received some metagenomic data from 16S rrna sequencing. The sequencing company claim to have removed primers, however not adapter sequences. Please note that the files have been demultiplexed.

In the platform, we have forward and reverse strand fastq files. For each fastq file, there is information regarding the adapter used, for each strand.

However, when I run the following command: `cat myFastQ_forward.fastq | grep 'adapter_forward'`, it returns nothing. My question is, shouldn't the adapters be identified successfully simply by looking for an exact match in the fastq files? If not, why?

Thanks,

# Answer
> 0 votes
It depends on the read layout they used, but typically when I hear anyone talk about removing adapters in paired-end sequencing, it's the occurrences of adapter sequences at the *ends* of each read from the *complementary strand's* adapter. For example see the figures on this IDT page about 16S sequencing. If a given sequence is short enough/number of sequencing cycles high enough, you'll end up seeing "past the ends" on each read.

For example, with a reverse read reverse-complemented and aligned with a forward read, you could maybe see something like this off the sequencer:

```
                       |-->
Forward:               BarcodefwdAmpliconBarcoderevAdapterrev...
Reverse:  ...AdapterfwdBarcodefwdAmpliconBarcoderev
                                               <--|

```

And then if the company already removed barcodes and primers at the start of the reads, maybe they gave you:

```
                                 |-->
Forward:                         AmpliconBarcoderevAdapterrev...
Reverse:  ...AdapterfwdBarcodefwdAmplicon
                                     <--|

```

So the material you'd want to trim is on the ends, and reverse-complemented. How much you see depends on the lengths involved so maybe you'll see partial adapters or none at all, and there may be sequencing error, too. See cutadapt's documentation about trimming 3' adapters and paired-end reads for more on that.

---
Tags: linux, metagenome, trimming, cutadapt, 16rrna
---


# Question
Title: Identifying Novel Species with ANI Taxonomy
I have sequenced a bacterial genome and after I checked the ANI scores, the scores against similar bacterias were 79 and lower. How can I confirm if this bacteria is a novel species/strain? Is there any threshold value for ANI score? I'm pretty sure if the ANI score is above 95, it indicates the same species.

Should I double-check this result by using other methods such as GGDC?

Thanks in advance!

# Answer
> 2 votes
## Short answer:

Yes.

## Long answer:

I suggest 4 methods for checking the genome.

1. Use minimap2 or a similar aligner to perform alignments with the most similar genomes, then generate a dot plot or other visualization of the regions that align to figure out what is going on. Here is one tool, I haven't used it, here is another. I generally use minimap2 to very quickly align, and then a quick plot in R using the PAF file.
2. Use k-mer hashing methods such as mash or sourmash to quickly compare to a large database. Large precomputed databases are available for download on the linked mash site. You can also use it for one-by-one comparisons.
3. Use QUAST to estimate more informative summary statistics for comparisons to the more similar genomes. QUAST has a rich set of possible visualizations that you can use. It also tries to identify the most closely related public genome using rDNA sequences, though that has been hit or miss for me.
4. Use a quality assessment tool such as CheckM or BUSCO to estimate the quality of your genome. THese also have utilities for taxonomic placement of your genome in case something bad happened (e.g. sample swap).

---
Tags: genome
---


# Question
Title: 10X low rate of correct barcodes was observed for the candidate chemistry choices for the input
I am testing a 10x fastq dataset , but cellranger count complains "An extremely low rate of correct barcodes was observed for all the candidate chemistry choices for the input", I have tried all the chemistry recommendation codes on the 10x genomics website, and I still get the same result, how can I solve this problem, thanks.

# Answer
I have solve this problem, cause the combined barcode is differernt from the common 10x barcode. when mapping the genome, need to add-- chemistry ARC-v1 when cellranger count.
> 2 votes

# Answer
Answer from @maximilian-press, converted from comments:

> I can't comment on 10X data specifically, however, I suggest writing a simple script or command that counts the frequency of each barcode, and reports back the N (10? 20? 100?) most frequent ones, which will allow you to investigate in more detail. This is what I have done in the case of suspected bardode swaps or reverse complement issues.
> 0 votes

---
Tags: single-cell, fastq
---


# Question
Title: Extracting Sequence Unique to a Certain Genome
I have two genomes with high similarity. But I found out that one genome has a longer sequence of about 200 kbp. I try to align the sequence with Mauve. How can I extract this unique 200 kbp?

# Answer
I've found the answer.

I've aligned these genomes (gbk files) with Mauve Contig Mover, and then I exported the Gap File into a csv file. Then I used python pandas to analyze the region that is unique to each strand from the csv file.
> 1 votes

---
Tags: sequence-alignment, genome, mauve
---


# Question
Title: Converting Control and Case vcf files into a format that can be processed by plink 1.9.0
I have two vcf files, one containing samples with TP53 mutations, and one containing samples with no TP53 mutations. The vcf without the mutations is the control and the vcf with the mutations is the case. How can I combine these two files into something that can be processed by the --assoc function in plink 1.9.0? Do I need to make them into one vcf? Or is it not necessary, and I can just use --assoc on each file after converting each individually?

# Answer
Both Plink 1.9 \[uses `--assoc`\] and Plink 2 \[uses `--glm`\] can process VCF files directly; no conversion needed.

I think that Plink will combine multiple files specified on the command line into one file, but if this is not the case then bcftools can be used to merge them.

As described in those input specifications, you can provide Plink with a `.fam` file (via `--fam`) to add in phenotype data.

Information about the `.fam` formats can be found here. It's a 6-column headerless space-separated format containing the following fields:

1. Family ID ('FID')
2. Within-family ID ('IID'; cannot be '0')
3. Within-family ID of father ('0' if father isn't in dataset)
4. Within-family ID of mother ('0' if mother isn't in dataset)
5. Sex code ('1' = male, '2' = female, '0' = unknown)
6. Phenotype value ('1' = control, '2' = case, '-9'/'0'/non-numeric = missing data if case/control)

The phenotype column can be quantitative as well; which is triggered when a numeric value is seen that is not one of the above qualitative codes.
> 1 votes

---
Tags: vcf, plink
---


title: Comparing homozygosity of k-mer plots

text:
# Question
Title: Comparing homozygosity of k-mer plots
Attached are two kmer plots from two closely related species. Is that safe to say that the one on the left has higher homozygosity than the one in the right k-mer plot, due to a low to almost flat blue peak (2x) in the left plot? My concern is that the blue peaks in the two plots are not in the same position: left plot at 4n, while right plot at 2n. I am aware of the grey area in the right plot indicating missing content in assembly, likely due to reads from multiple individuals (correct me if I am wrong).

> **Model for kmer plot on the left**
>
> ```
> Formula: y_transform ~ x^transform_exp * length * predict2_0(r1, k, d,
>     kmercov, bias, x)
>
> Parameters:
>          Estimate Std. Error t value Pr(>|t|)
> d       4.019e-02  1.262e-03   31.84   <2e-16 ***
> r1      1.014e-02  1.664e-04   60.94   <2e-16 ***
> kmercov 1.996e+01  2.452e-02  814.12   <2e-16 ***
> bias    1.900e+00  2.495e-02   76.17   <2e-16 ***
> length  4.553e+08  8.993e+05  506.34   <2e-16 ***
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 4024000 on 1995 degrees of freedom
>
> Number of iterations to convergence: 8
> Achieved convergence tolerance: 1.49e-08
>
> ```

> **Model for kmer plot on the right**
>
> ```
> Formula: y_transform ~ x^transform_exp * length * predict2_0(r1, k, d,
>     kmercov, bias, x)
>
> Parameters:
>          Estimate Std. Error t value Pr(>|t|)
> d       4.399e-02  1.407e-03   31.27   <2e-16 ***
> r1      1.519e-02  1.002e-04  151.64   <2e-16 ***
> kmercov 2.069e+01  2.416e-02  856.53   <2e-16 ***
> bias    1.697e+00  2.099e-02   80.88   <2e-16 ***
> length  5.069e+08  1.056e+06  480.10   <2e-16 ***
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 4878000 on 1995 degrees of freedom
>
> Number of iterations to convergence: 7
> Achieved convergence tolerance: 1.49e-08
>
> ```

Both the GenomeScope estimated genome sizes are smaller than the assembly (both have microorganisms filtered). The genome size of the left species is ~100Mb smaller than the one on the right, consistently in GenomeScope and the assembly.

## GenomeScope Plots:

# Answer
I would say yes, but it is not that easy to see it on a histogram decomposed to overlapping distributions. While the decomposition to different ditributions is assembly dependent, the overall k-mer histogram is actually directly derived from the sequencing data.

And these k-mer histograms can be used to estimate heterozygosity. For example, using genomescope. You can install the tool locally, or fit the model on a webserver (which seems to be down at the moment, but it should be back on at some point). For more details check the paper. If you get the tool running, you can estimate heterozygosity in the two samples to see how much more/less heterozygous they are.

P.S. I am always thinking about heterozygosity, which is supposed an equivalent of 1 - homozygosity, so I hope the answer will make sense to you.

--- edit ---

so, looking at the gneomescope models and knowing that genome size estimate matches the assembly size, I think it's fair to conclude the right genome is more heterozygous than the left one. However, keep in mind that it is not super great sequencing dataset, as the error peak and heterozygous peak are not clearly separated, hence the estimates of genomic features might be less accurate.
> 1 votes

---
Tags: genome, assembly, k-mer
---


# Question
Title: How to find an unidentified splicing variant of a protein?
I have some data indicating there might be a splicing variant of the Arabidopsis Thaliana protein I'm studying that has not been identified. Is there a database of i.e. RNA sequences (transcriptome?) that I could search for it, knowing a part of its sequence?

# Answer
> 1 votes
There's not going to be a specific database that's particularly thorough for this. The best you can do is one of the following:

1. Search through TAIR (https://www.arabidopsis.org/servlets/Search?action=new\_search&type=protein)
2. Blast SRA with the full length sequence ( https://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastn&BLAST\_PROGRAMS=megaBlast&PAGE\_TYPE=BlastSearch&BLAST\_SPEC=SRA&SHOW\_DEFAULTS=on). This will unfortunately pull up a LOT of related hits for every other related isoform.
3. Download some public RNA-seq data from SRA/ENA and see if it has supporting data.

You're probably best off doing all of these, where you use blast to find relevant datasets that you then align. My guess is that TAIR won't have everything, since this has probably been seen before but not reported since the experiment wasn't intending to find new proteins.

---
Tags: rna-seq, transcriptome, rna-splicing
---


# Question
Title: Length of Contigs in Transcriptome and Whole Genome Assembly
Why are there shorter contigs from transcriptome assembly than from a whole genome assembly? I know the difference between transcriptome and genome, but don't really understand what contigs are in the context of sequencing in bioinformatics, and why are the contigs shorter in transcriptome assembly than in whole genome assembly. Anyone mind explaining? Much appreciated.

# Answer
The longest possible transcriptome contig reflects the longest possible transcript. In the human genome, that's possibly Titin, which is ~35K aa (that's ~105Kb) and will be longer pre-splicing.

The longest possible genome contig reflects the longest chromosome. In the human genome that's chr1 at ~250Mb.

In practise your contigs/scaffolds will be a lot shorter for various reasons (e.g. hard to sequence the telomeres and repeats) but it illustrates the difference since all the chromosomes consist of contiguous DNA but many separate RNA transcripts (and not all DNA is transcribed).
> 3 votes

---
Tags: assembly, transcriptome, genome-sequencing
---


# Question
Title: List of all Latin word parts used in Biological Taxonomy?
I seem to have determined that there are 287,614 unique Latin taxonomic names used in the biological taxonomy according to the Integrated Taxonomic Information System (ITIS) taxonomy taxon.txt file. First of all, that is an enormous number of unique names, I would have hoped they could be reused to a greater degree but I guess not. There are 90k some odd genera names, the rest are species names, or certain names might be duplicated throughout the complete taxonomy.

What I would like to do is take the taxonomy Latin names and convert them into their English equivalents. Not their "common names", as that can't be done by looking at the Latin name alone. But by taking their Latin name and making it into the equivalent English literal meaning. Wikipedia has a nice list on commonly used taxonomic affixes, and another on Latin and Greek words commonly used in systematic names, however a quick glance at the list and they both are missing `-bacter-` as one infix, for a single example, or `vorans`, or many many others. The Wikipedia list seems to only contain a small fraction of the names used in the actual complete taxonomy.

What I'm wondering is if there exist anywhere online in text format (so it can be copy-pasted) a complete list of all the Latin affixes (prefixes, infixes, suffixes, etc.) so I can accomplish my goal. The goal being, programmatically convert the Latin names to their literal English equivalent. I would split apart the Latin name based on some database of Latin affixes, and then plop out the English word matching the affix. It should take only a few seconds for a computer to do this transformation, assuming the list mapping Latin names to English names exists already.

Does anything like this exist anywhere on the web?

This calls itself the "worlds largest affix directory", yet it still doesn't have these common ones... Here's another Wikipedia page, perhaps this is as good as it is going to get?

# Answer
Maybe you can use NCBI Taxonomy instead, which is a full hierarchical taxonomy with \>700K taxa. You can use Entrez utilities to programmatically query NCBI taxonomy to retrieve Latin or common names as well as other metadata.
> 1 votes

---
Tags: database, taxonomy
---


# Question
Title: How to retrieve fasta sequence after local blast?
I have created a Blast database using a reference genome. Then, I have performed a local blast search in command line using a gene of interest. I have obtained some hits with the usual Blasting information. Now, I want to extract the exact matching sequence from the blast search with their corresponding start and stop position. I have a `.gff` file for the genes as well for the reference genome. Can you give me any suggestions on how to do that? Thanks a lot.

# Answer
FWIW, In the Sequenceserver BLAST software, you can just click on "download alignment" to extract the aligning region.
> 0 votes

# Answer
My usual approach for extracting sequence from a reference genome is to use `samtools faidx`:

```
# create genome index; slow; only do once per reference genome
samtools faidx reference_genome.fa
# extract sequence for a specific match
samtools faidx reference_genome.fa chromosome_name1:posStart1-posEnd1 > match1.fa
samtools faidx reference_genome.fa chromosome_name2:posStart2-posEnd2 > match2.fa

```

You haven't mentioned how many BLAST results you have, or what format they're in, so I'm not sure if it's practical to do this manually, or if it needs to be automated in some way.
> 1 votes

---
Tags: gene, ngs, blast
---


title: maxEE parameter from DADA2 package in R

text:
# Question
Title: maxEE parameter from DADA2 package in R
I have just started the QC steps from the dada2 pipeline, and have failed to find a detailed explanation of what the maxEE argument entails.

I have surfed many forums, as well as the details given by the creators of the package, but they are lacking in detail. All it says is that:

```
After truncation, reads with higher than maxEE "expected errors" will be discarded. Expected
errors are calculated from the nominal definition of the quality score: EE = sum(10^(-Q/10))

```

What is Q in this instance? Is it the Quality score obtained from the `plotQualityProfile` function?

In the tutorial, it states that:

```
The standard filtering parameters are starting points, not set in stone. If you want to speed up
downstream computation, consider tightening maxEE. If too few reads are passing the filter,
consider relaxing maxEE, perhaps especially on the reverse reads (eg. MaxEE = c (2,5)), and
reducing the truncLen to remove low quality tails.

```

What is 2, and 5 in this instance? What does an expected error of 2, or 5, actually mean?

Thanks

# Answer
> 1 votes
The previous paragraph in the tutorial provides a link to the paper.

> We’ll use standard filtering parameters: maxN=0 (DADA2 requires no Ns), truncQ=2, rm.phix=TRUE and maxEE=2. The maxEE parameter sets the maximum number of “expected errors” allowed in a read, which is a better filter than simply averaging quality scores.

---
Tags: r, microbiome, 16rrna
---


title: What are Dummy fastqs?

text:
# Question
Title: What are Dummy fastqs?
I had to do a file transfer for in our system to run from fastqs and I came across the term "dummy fastq", not sure what exactly it means and what is the purpose of it in the workflow. Can someone answer?

# Answer
Probably some made-up data solely for the sake of testing a workflow but actually impossible to know what random files are sitting on your servers.
> 2 votes

---
Tags: ngs, linux, cancer
---


title: How to convert gff to gtf?

text:
# Question
Title: How to convert gff to gtf?
My annotation file is in .gff format. I would like to convert it to .gtf format or to know if there is a way to directly download the annotation file in .gtf format?

I am working on sequences from the P Falicparum 3D7.

# Answer
I would suggest to use `agat_convert_sp_gff2gtf.pl` from AGAT because you loose information with `gffread`.
e.g here a gff example:

```
##gff-version 3
scaffold625 maker   gene    337818  343277  .   +   .   ID=CLUHARG00000005458;Name=TUBB3_2
scaffold625 maker   transcript  337818  343277  .   +   .   ID=CLUHART00000008717;Parent=CLUHARG00000005458
scaffold625 maker   CDS 337915  337971  .   +   0   ID=CLUHART00000008717:cds;Parent=CLUHART00000008717
scaffold625 maker   CDS 340733  340841  .   +   0   ID=CLUHART00000008717:cds;Parent=CLUHART00000008717
scaffold625 maker   CDS 341518  341628  .   +   2   ID=CLUHART00000008717:cds;Parent=CLUHART00000008717
scaffold625 maker   CDS 341964  343033  .   +   2   ID=CLUHART00000008717:cds;Parent=CLUHART00000008717
scaffold625 maker   exon    337818  337971  .   +   0   ID=CLUHART00000008717:exon1;Parent=CLUHART00000008717
scaffold625 maker   exon    340733  340841  .   +   0   ID=CLUHART00000008717:exon2;Parent=CLUHART00000008717
scaffold625 maker   exon    341518  341628  .   +   2   ID=CLUHART00000008717:exon3;Parent=CLUHART00000008717
scaffold625 maker   exon    341964  343277  .   +   2   ID=CLUHART00000008717:exon4;Parent=CLUHART00000008717
scaffold625 maker   five_prime_UTR  337818  337914  .   +   .   ID=CLUHART00000008717:five_prime_utr;Parent=CLUHART00000008717
scaffold625 maker   three_prime_UTR 343034  343277  .   +   .   ID=CLUHART00000008717:three_prime_utr;Parent=CLUHART00000008717

```

with gffread you get:

```
scaffold625 maker   transcript  337818  343277  .   +   .   transcript_id "CLUHART00000008717"; gene_id "CLUHARG00000005458";
scaffold625 maker   exon    337818  337971  .   +   .   transcript_id "CLUHART00000008717"; gene_id "CLUHARG00000005458";
scaffold625 maker   exon    340733  340841  .   +   .   transcript_id "CLUHART00000008717"; gene_id "CLUHARG00000005458";
scaffold625 maker   exon    341518  341628  .   +   .   transcript_id "CLUHART00000008717"; gene_id "CLUHARG00000005458";
scaffold625 maker   exon    341964  343277  .   +   .   transcript_id "CLUHART00000008717"; gene_id "CLUHARG00000005458";
scaffold625 maker   CDS 337915  337971  .   +   0   transcript_id "CLUHART00000008717"; gene_id "CLUHARG00000005458";
scaffold625 maker   CDS 340733  340841  .   +   0   transcript_id "CLUHART00000008717"; gene_id "CLUHARG00000005458";
scaffold625 maker   CDS 341518  341628  .   +   2   transcript_id "CLUHART00000008717"; gene_id "CLUHARG00000005458";
scaffold625 maker   CDS 341964  343033  .   +   2   transcript_id "CLUHART00000008717"; gene_id "CLUHARG00000005458";

```

while with AGAT you get

```
scaffold625 maker   gene    337818  343277  .   +   .   ID CLUHARG00000005458 ; Name TUBB3_2 ; gene_id CLUHARG00000005458
scaffold625 maker   mRNA    337818  343277  .   +   .   ID CLUHART00000008717 ; Parent CLUHARG00000005458 ; gene_id CLUHARG00000005458 ; transcript_id CLUHART00000008717
scaffold625 maker   exon    337818  337971  .   +   0   ID "CLUHART00000008717:exon1"  ; Parent CLUHART00000008717 ; gene_id CLUHARG00000005458 ; transcript_id CLUHART00000008717
scaffold625 maker   exon    340733  340841  .   +   0   ID "CLUHART00000008717:exon2"  ; Parent CLUHART00000008717 ; gene_id CLUHARG00000005458 ; transcript_id CLUHART00000008717
scaffold625 maker   exon    341518  341628  .   +   2   ID "CLUHART00000008717:exon3"  ; Parent CLUHART00000008717 ; gene_id CLUHARG00000005458 ; transcript_id CLUHART00000008717
scaffold625 maker   exon    341964  343277  .   +   2   ID "CLUHART00000008717:exon4"  ; Parent CLUHART00000008717 ; gene_id CLUHARG00000005458 ; transcript_id CLUHART00000008717
scaffold625 maker   CDS 337915  337971  .   +   0   ID "CLUHART00000008717:cds"  ; Parent CLUHART00000008717 ; gene_id CLUHARG00000005458 ; transcript_id CLUHART00000008717
scaffold625 maker   CDS 340733  340841  .   +   0   ID "CLUHART00000008717:cds"  ; Parent CLUHART00000008717 ; gene_id CLUHARG00000005458 ; transcript_id CLUHART00000008717
scaffold625 maker   CDS 341518  341628  .   +   2   ID "CLUHART00000008717:cds"  ; Parent CLUHART00000008717 ; gene_id CLUHARG00000005458 ; transcript_id CLUHART00000008717
scaffold625 maker   CDS 341964  343033  .   +   2   ID "CLUHART00000008717:cds"  ; Parent CLUHART00000008717 ; gene_id CLUHARG00000005458 ; transcript_id CLUHART00000008717
scaffold625 maker   five_prime_UTR  337818  337914  .   +   .   ID "CLUHART00000008717:five_prime_utr"  ; Parent CLUHART00000008717 ; gene_id CLUHARG00000005458 ; transcript_id CLUHART00000008717
scaffold625 maker   three_prime_UTR 343034  343277  .   +   .   ID "CLUHART00000008717:three_prime_utr"  ; Parent CLUHART00000008717 ; gene_id CLUHARG00000005458 ; transcript_id CLUHART00000008717

```

So you can see as example that `Name=TUBB3_2` from the gene feature disappeared with `gffread`.
> 6 votes

# Answer
You can use `gffread` to convert `gff` to `gtf2`, below is from the manual:

> In order to see the GTF2 version of the same transcripts, the -T option should be added:
>
> ```
> gffread -E annotation.gff -T -o- | more
>
> ```
>
> The examples above also show that gffread can be used to convert a file between GTF2 and GFF3 file formats.
> 5 votes

# Answer
Actually I started a mini-review about the tools to do such conversion.
You can find it here: https://github.com/NBISweden/GAAS/blob/master/annotation/knowledge/gff\_to\_gtf.md

As I mentioned earlier they do not behave all the same way...
> 1 votes

# Answer
You can try bioinfokit (https://github.com/reneshbedre/bioinfokit) in python.

from bioinfokit.analys import gff

gff.gff_to_gtf(file="yourfile.gff3")
# The converted gtf file will be saved in the same working directory.
```

> 0 votes

---
Tags: rna-seq, gtf, gff
---


# Question
Title: Using different "principal components" in UMAP visualizations
In certain UMAP visualization, it appears that the visualizations show the first UMAP "principal component" plotted against the second. For example, like this:

Apparently, this lets us use the first and second UMAP "principal components" to tell apart the different cell clusterings in a dataset. I'm curious about if we were to use e.g. the second and the *third* UMAP "principal components" instead - would this let us tell apart different cell clusterings according to the second and third features that explain those clusterings' largest sources of variation? If someone could explain this to me, I would be very grateful, and if you would like to correct whatever misunderstandings I have about UMAP as well, that would be great too - I'd be the first one to tell you that I don't understand it well!

# Answer
You can plot UMAP3 too in 3 dimensions, but that's about the limit for viz. You could plot UMAP2 vs. UMAP3, UMAP3 vs. UMAP4 etc. but I'm not sure how informative it is. Instead, you can focus on other analyses like clustering, differential expression, cell type annotation, marker gene analysis, and other machine learning approaches etc. Ultimately, UMAP/t-SNE are nice tools for viz but they can be over interpreted- it is inherently impossible to summarise very highly dimensional data in just a few dimensions without losing information- see here for an (overly) provocative take.
> 1 votes

---
Tags: scrnaseq, seurat, visualization, umap
---


# Question
Title: Remove variable sequence component within a tree text file
I have a gene tree file of 436 orthologue genes from 6 species. I want to remove unwanted extensions as it looks massy after visualization. My file looks like:

(TRINITY\_Clupea\_DN5452\_c0\_g1\_i1.p1:0.0824467436,TRINITY\_Engraulis\_DN43599\_c0\_g1\_i1.p1:0.1634781085)100:0.0876433106,TRINITY\_Sardina\_DN15766\_c0\_g1\_i2.p1:0.0164132018)..................

What i need:

(Clupea\_DN5452:0.0824467436,Engraulis\_DN43599:0.1634781085)100:0.0876433106,Sardina\_DN15766:0.0164132018)..................

As "TRINITY" is identicial, i can remove it using sed. But after the species name ids are not identical. And i only need the 2nd and 3rd part of the identifier.

It will be helpful for me if i get some suggestion. Thanks.

# Answer
Lots of ways to do this, but here's one way using :

```
sed 's/TRINITY_\([^_]\+\)_\([^_]\+\)_[^:]\+/\1_\2/g' file.tree

```

The above applies a substitution to capture the first two fields using parentheses (which need to be escaped). Each field captures any character not an underscore one or more times (using the "+" symbol, which also needs to be escaped). Finally, we also match any character not a colon one or more times. The replacement string can be just the first two capture groups separated with an underscore.
> 1 votes

# Answer
If you want to remove component/gene/isoform parts from a Trinity ID, you can do `s/_c[0-9]+_.*?:/:/` to remove anything between the first c with only numbers after it up to the next colon. Combining this with the TRINITY removal at the start, it'd be something like this:

```
sed 's/TRINITY_\(.*?\)_c[0-9]+_.*?:/\1:/g' file.tree

```
> 0 votes

---
Tags: phylogenetics, identifiers, trinity
---


title: Why does randomized motif search work?

text:
# Question
Title: Why does randomized motif search work?
I'm looking for intuition for why a randomized motif search works. My current thinking is as follows:

We are selecting many random kmers from our DNA sequences. The chosen kmers will bias the profile matrix to selecting kmers like them.

Given any particular k-mer chosen, there are two possibilities:

1. We've selected a meaningless kmer, so the profile will bias selecting random (ie high entropy) kmers and the score will stop improving soon.
2. We've captured a portion of the true motif. Since the motif isn't random (ie doesn't have high entropy) subsequent iterations using this profile matrix will pick motifs closer to the real motif and the score will keep improving.

Since we give this algorithm many tries and at some point it will eventually arrive at possibility 2 and come to a good conclusion.

Is this a correct understanding of why randomized motif search works? Is there a more rigorous way to prove this?

# Answer
> 2 votes
## Update

With regard to comment, the proof of effectiveness of Monte Carlo algorithms like this is a little thorny- if you look at the wikipedia page, you will see that a proof of statistical consistency for any Monte Carlo algorithms took a long time:

> From 1950 to 1996, all the publications on Sequential Monte Carlo methodologies, including the pruning and resample Monte Carlo methods introduced in computational physics and molecular chemistry, present natural and heuristic-like algorithms applied to different situations without a single proof of their consistency, nor a discussion on the bias of the estimates and on genealogical and ancestral tree based algorithms. The mathematical foundations and the first rigorous analysis of these particle algorithms were written by Pierre Del Moral in 1996.

Here is a possibly more general reference by Del Moral in 1998. The paper is somewhat dense and I am not an expert in the field so I don't think I'll do it justice, but I think that we could probably figure out a similar mathematical proof for random motif search given time and the right problem formulation.

I sense that your concern might be with Monte Carlo algorithms generally rather than random motif search specifically, so hopefully that gets close to your question.

It does look like people have done heuristic analysis of similar algorithms, but they are not formal proofs.

## Original answer

I believe that you are more or less correct.

For example, this presentation's walkthrough of the algorithm (slides 35-36) specifically refers to **Greedy** randomized profile motif searches. By making the greediness of the algorithm explicit, we can see that there is no claim of finding an optimal solution in any particular iteration, thus the need for multiple retries to get somewhere close to approximate solution.

One additional point made in that presentation is that this algorithm works best when coupled to approaches such as Gibbs sampling for targeting the sampled k-mer chain in regions with higher scores (less entropic k-mers).

A different learning resource with similar content (and a little more intuition) is available here.

---
Tags: sequence-analysis, k-mer, algorithms, motifs
---


# Question
Title: Different results in differential expression analysis of microarray data
I am performing differential gene expression analysis to microarray data for type 2 diabetes donors and nondiabetic donors. When I run the code I get some different results in each time (about 50 or so different genes in each run). Is this normal? I am submitting a project and it could be very problematic if my results change in each run, since I will perform downstream analysis for the DE genes and use web tools such as DAVID for further analysis.

Is there a way to keep the results from changing in each run?

```
g <- getGEO("GSE20966",GSEMatrix = TRUE)
g<- g[[1]]
exprs(g)
exprs(g) <- log2(exprs(g))

#Diffrential gene exp

library(limma)
design <- model.matrix(~0+sampleInfo$characteristics_ch1.1)

colnames(design) <- c("NonDiabetic","T2D")

summary(exprs(g))

exp <- na.omit(exprs(g))
cutoff <- median(exp)

## TRUE or FALSE for whether each gene is "expressed" in each sample
is_expressed <- exprs(g) > cutoff

## Identify genes expressed in more than 2 samples

keep <- rowSums(is_expressed) > 2

## check how many genes are removed / retained.
table(keep)
keep<- na.omit(keep)

## subset to just those expressed genes
g<- g[keep,]
fit <- lmFit(exp, design)

head(fit$coefficients)
contrasts <- makeContrasts(T2D - NonDiabetic, levels=design)
fit2 <- contrasts.fit(fit, contrasts)
fit2 <- eBayes(fit2)
topTable(fit2)
topTable(fit2, coef=1)
decideTests(fit2)
table(decideTests(fit2))
aw <- arrayWeights(exp,design)
aw
fit <- lmFit(exp, design,
            weights = aw)
contrasts <- makeContrasts(T2D - NonDiabetic, levels=design)
fit2 <- contrasts.fit(fit, contrasts)
fit2 <- eBayes(fit2)
g <- na.omit(g)
anno <- fData(g)
anno
anno <- select(anno,`Gene Symbol`,ENTREZ_GENE_ID)
fit2$genes <- anno
topTable(fit2)
full_results <- topTable(fit2, number=Inf)
full_results <- tibble::rownames_to_column(full_results,"ID")
library(ggplot2)
dev.off()
ggplot(full_results,aes(x = logFC, y=B)) + geom_point()
p_cutoff <- 0.05
fc_cutoff <- 1

full_results %>%
  mutate(Significant = adj.P.Val < p_cutoff, abs(logFC) > fc_cutoff )
 %>%
  mutate(Rank = 1:n(), Label = ifelse(Rank < topN, anno$`Gene
 Symbol`,"")) %>%
  ggplot(aes(x = logFC, y = B, col=Significant,label=Label)) +
geom_point() + geom_text_repel(col="black")

 p_cutoff <- 0.05
 fc_cutoff <- 1

filter(full_results, adj.P.Val < 0.05, abs(logFC) > 1)
 library(pheatmap)
topN <- 20

top<-150
gene_names <- mutate(full_results, Rank = 1:n()) %>%
  filter(Rank < top) %>%
  pull(Gene.Symbol)

```

# Answer
DE results do not change from one run to another. The code given in your question will give identical results each time you run it.
> 2 votes

---
Tags: r, bioconductor, differential-expression, gene-expression, microarray
---


# Question
Title: Determining what RNAseq data is filtered on volcano plot
I am using RNA seq data to analyze genes via a volcano plot comparing differential gene expression of bacteria with and without antibiotic in R. After having created my plot, I am unsure why some of my values which read "0" TPM in one condition and a TPM value that's not "0" in the other condition were not determined to be differentially expressed. Some of the genes on my volcano plot have this difference in TPM and show up as significant on my plot while others with this "0" value difference are not considered to be significant according to my plot.

**Sample of my data:** (UO1\_D712\_##### represents the locus number, the top column represents the different replicates where "1" is untreated and "2" is treated with antibiotic, the letters next to the numbers represent a different technical replicate):

```
structure(list(`1` = c("U01_D712_00001",

"U01_D712_00002", "U01_D712_00003",
"U01_D712_00004", "U01_D712_00005", "U01_D712_00006", "U01_D712_00007",
"U01_D712_00008", "U01_D712_00009", "U01_D712_00010", "U01_D712_00011",
"U01_D712_00012", "U01_D712_00013", "U01_D712_00014", "U01_D712_00015",
"U01_D712_00016", "U01_D712_00017", "U01_D712_00018", "U01_D712_00019",
"U01_D712_00020", "U01_D712_00021", "U01_D712_00022", "U01_D712_00023",
"U01_D712_00024", "U01_D712_00025", "U01_D712_00026", "U01_D712_00027",
"U01_D712_00028", "U01_D712_00029", "U01_D712_00030"), `1a` = c("6.590456502",
"6.741694688", "7.110342585", "6.299482433", "2.77173982", "2.470330508",
"3.827125186", "6.267229842", "5.524708663", "2.657913228", "2.87209859",
"7.479820548", "4.980185572", "4.210955199", "0", "4.492492822",
"3.611091371", "3.714270433", "7.455036914", "7.045203025", "6.061860857",
"2.925268313", "6.544077039", "5.747318013", "9.97083444", "9.22523089",
"4.20205383", "6.097040679", "2.621192351", "0"), `1b` = c("6.544454427",
"6.601489488", "7.134224619", "5.814043553", "3.280958379", "2.649180803",
"3.860542083", "6.256648363", "5.380427766", "2.581027705", "3.016132165",
"7.405329447", "4.701503289", "4.073814818", "0", "4.304196924",
"3.515977329", "3.843535649", "7.342972625", "6.966606769", "6.122878624",
"3.007522306", "6.495797641", "5.965431621", "9.828050269", "9.219563915",
"4.065778989", "6.105331066", "3.061209408", "0"), `1c` = c("6.608196006",
"6.743010138", "7.102600793", "6.146518601", "3.555184202", "2.364971542",
"4.034053983", "6.158523627", "5.656051812", "2.658660735", "3.054717455",
"7.392164473", "4.950953264", "4.277770477", "0", "4.507936666",
"3.794842979", "3.610794578", "7.471646548", "7.104792624", "6.484767016",
"3.071205184", "6.584425715", "6.20466015", "9.986342122", "9.282943758",
"4.179958213", "6.219551653", "2.984738345", "0"), `1d` = c("6.547155382",
"6.558328892", "6.992501615", "6.449558793", "3.059801464", "2.418800257",
"3.96498952", "6.013208538", "5.279919645", "2.893295085", "1.750510471",
"7.408735671", "4.9425624", "3.804549986", "0", "4.174835979",
"3.806006888", "3.570390524", "7.641137006", "6.976672494", "6.363030106",
"3.083061726", "6.300910093", "6.007490342", "9.926316442", "9.09671588",
"4.320556917", "6.153860107", "2.877230446", "0"), `1e` = c("6.626724417",
"6.577345176", "7.156821278", "6.296873411", "3.618089702", "2.394444986",
"4.129376392", "6.011715246", "5.31197869", "2.00754706", "2.695493528",
"7.538910448", "4.606060035", "3.909472643", "0", "4.346616047",
"3.468681284", "3.338231445", "7.559599613", "7.1527452", "6.232923513",
"3.108624209", "6.535435309", "6.12922864", "10.13108497", "9.310331313",
"3.959568571", "6.182335537", "2.902736258", "0"), `1f` = c("6.419219179",
"6.650459302", "7.319125725", "5.570869357", "2.962845933", "2.55903176",
"4.087597573", "5.995610538", "5.386268651", "2.750800859", "2.416572678",
"7.579148955", "3.952633067", "3.615227674", "0", "4.562838935",
"3.76942104", "3.920096905", "7.935320749", "7.501470652", "6.13700099",
"3.123910608", "6.035971952", "5.706235015", "9.254751395", "8.379630979",
"4.51391973", "5.6890651", "2.43285316", "0"), `1g` = c("6.553221221",
"6.633949852", "7.182305386", "5.769365973", "2.721354972", "1.668390466",
"4.148367057", "6.240883382", "5.458877133", "1.733842637", "2.723803355",
"7.522249899", "4.149567197", "3.780763096", "0", "4.496306813",
"3.645643535", "3.851001768", "7.678552875", "7.283411279", "6.591585956",
"2.879345378", "6.389427003", "5.911222165", "9.851084493", "9.084575304",
"4.272587776", "5.974762147", "2.98852705", "0"), `2a` = c("9.769887737",
"4.550226652", "1.869021464", "4.944848987", "7.9678549", "2.682865013",
"8.495575559", "2.234521659", "3.667196316", "9.180445037", "3.210107621",
"7.21523691", "5.714579923", "5.423986751", "9.118981459", "8.635701597",
"2.742889473", "3.712618983", "8.006057144", "4.999541279", "10.54351774",
"5.880978085", "7.145433526", "6.982416661", "9.339651188", "5.360835327",
"4.699680905", "3.423826225", "6.408271885", "5.038170992"),
    `2b` = c("10.26397519", "1.945664005", "2.086763158", "2.800904763",
    "8.583418657", "1.536094563", "9.32057547", "3.10685839",
    "4.224502319", "10.1252842", "1.811175407", "5.714439316",
    "6.039142559", "3.833361174", "9.757360286", "9.565906731",
    "1.523640473", "2.315033488", "6.312524363", "4.889986456",
    "10.23020108", "4.848727685", "7.533256999", "7.138160378",
    "10.30380331", "3.955469283", "3.167940742", "3.599655687",
    "4.828945262", "3.701043054"), `2c` = c("9.478481216", "2.789289131",
    "3.393949527", "3.754810933", "8.710777154", "1.806170784",
    "9.150005253", "2.612275457", "4.961073313", "9.802701699",
    "2.933183115", "6.532384958", "6.919449225", "4.432699799",
    "9.715063475", "9.265691356", "1.412064593", "3.330131873",
    "6.665979896", "5.158526421", "9.417365584", "4.899531204",
    "8.173459354", "7.271400938", "9.813068613", "4.384622077",
    "3.700645365", "4.457874829", "5.649440022", "3.531010379"
    ), `2d` = c("8.199795497", "2.565711524", "2.202287889",
    "3.856354444", "7.380224849", "2.192476466", "8.14446837",
    "1.144258862", "3.31447122", "8.713146629", "2.697890381",
    "6.304428859", "5.745291803", "4.898396114", "9.173362747",
    "8.339933849", "3.159678152", "4.094587234", "7.608649692",
    "5.280206424", "10.34630403", "5.098585806", "7.262400625",
    "6.150190905", "9.316698845", "5.073027993", "4.695003229",
    "2.485847024", "5.545300465", "4.350571411"), `2e` = c("8.74935033",
    "3.739489484", "0.217205045", "4.413657999", "8.745525588",
    "3.657060086", "9.279834921", "2.898951179", "4.282874018",
    "9.610485827", "3.561102455", "7.228334332", "6.388491443",
    "4.389908652", "8.781086564", "9.178866581", "2.374603596",
    "3.961037408", "7.864369809", "3.654728044", "10.15284858",
    "5.894439123", "7.68020282", "7.12243523", "9.998637438",
    "5.092174395", "4.111530392", "3.776835632", "5.624523213",
    "4.095011377"), `2f` = c("7.648310926", "4.345557215", "1.986576876",
    "4.99426288", "7.087937177", "2.810917253", "7.77706637",
    "2.62822773", "3.581811188", "8.470225989", "3.335757437",
    "7.416094847", "5.208841128", "5.536128034", "8.255571138",
    "7.993319997", "1.9209089", "3.573861828", "7.318814519",
    "5.233804806", "11.05855833", "6.247720809", "6.673407583",
    "6.029960625", "8.806867591", "5.459208493", "4.001428729",
    "3.609936979", "5.876522973", "4.652839671"), `2g` = c("7.555235468",
    "3.892899549", "1.726443458", "5.304546796", "7.039588042",
    "3.027235295", "7.703852207", "1.753183519", "2.909288568",
    "8.385169315", "3.902707541", "7.523315081", "4.978364017",
    "5.49103181", "8.096218606", "7.944822989", "2.352609608",
    "4.155433517", "7.227355741", "5.532668321", "11.24946953",
    "6.159185473", "6.443203375", "5.931761874", "8.7421732",
    "5.502000205", "4.652883503", "3.458323017", "6.566487449",
    "4.89790353")), row.names = c(NA, 30L), class = "data.frame")

```

**Design matrix:**

```
structure(list(SampleID = c("1a ", "1b", "1c",

"1d", "1e", "1f",
"1g", "2a", "2b", "2c", "2d", "2e", "2f", "2g"), Group = c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), Replicate = c("a",
"b", "c", "d", "e", "f", "g", "a", "b", "c", "d", "e", "f", "g"
)), row.names = c(NA, 14L), class = "data.frame")

```

Here is how I wrote my code to sample the RNA seq data (note ~/Documents/VOLCANO/R10LB\_0ugvs0.5ug.csv refers to the the sample of data above and ~/Documents/VOLCANO/R10LB\_0vs0.5\_designmatrix.csv refers to the above design matrix):

```
R10LB_0vs0.5 <- read.csv2("~/Documents/VOLCANO/R10LB_0ugvs0.5ug.csv", sep=",", check.names = F)
R10LB_0vs0.5 <- janitor::remove_empty(R10LB_0vs0.5, which = "cols")
R10LB_0vs0.5_design <- read.csv2("~/Documents/VOLCANO/R10LB_0vs0.5_designmatrix.csv", sep=",")

rownames(R10LB_0vs0.5) <- R10LB_0vs0.5$`1`
R10LB_0vs0.5 <- R10LB_0vs0.5[,-1]

#colnames(R10LB_0ugvs0.5ug) <- R10LB_0ugvs0.5ug[1,]
#R10LB_0ugvs0.5ug <- R10LB_0ugvs0.5ug[-1,]

R10LB_0vs0.5 <- data.matrix(R10LB_0vs0.5)

R10LB_0vs0.5 <- DGEList(counts = R10LB_0vs0.5, group = R10LB_0vs0.5_design$Group)
str(R10LB_0vs0.5)

Group <- as.vector(as.character(R10LB_0vs0.5_design$Group))
Replicate <- as.vector((R10LB_0vs0.5_design$Replicate))
R10LB_0vs0.5_designmatrix <- model.matrix(~0+Group+Replicate)

keep <- filterByExpr(R10LB_0vs0.5)
R10LB_0vs0.5_filter <- R10LB_0vs0.5[keep, keep.lib.sizes = FALSE]

R10LB_0vs0.5_Disp <- estimateDisp(R10LB_0vs0.5_filter, R10LB_0vs0.5_designmatrix)

CONTRASTS <- makeContrasts( Group1vs2 = Group1-Group2,
                            levels = R10LB_0vs0.5_designmatrix)

R10LB_0vs0.5fit <- glmFit(R10LB_0vs0.5_Disp, contrast = CONTRASTS[,1])
R10LB_0vs0.5lrt <- glmLRT(R10LB_0vs0.5fit)
R10LB_0vs0.5TT <- topTags(R10LB_0vs0.5lrt, n=nrow(R10LB_0vs0.5_Disp))
write.csv(R10LB_0vs0.5TT, file = "R10LB_0vs0.5_comparison.csv")
saveRDS(R10LB_0vs0.5TT, file = "R10LB_0vs0.5_comparison.RDS")

```

**Volcano plot:**

```
ggplot(data = R10LB_0vs0.5_comparison, aes(x = logFC, y = PValue)) + geom_point()

p4 <- ggplot(data = R10LB_0vs0.5_comparison, aes(x = logFC, y = -log10(PValue))) + geom_point() + theme_minimal()

#Add vertical lines for logFC thresholds, and one horizontal line for the p-value threshold
p5 <- p4 + geom_vline(xintercept = c(-0.6, 0.6), col = "red") + geom_hline(yintercept = -log10(0.05), col = "red")

# add a column of NAs
R10LB_0vs0.5_comparison$diffexpressed <- "NO"
# if logFC > 0.6 and PValue < 0.05, set as "UP"
R10LB_0vs0.5_comparison$diffexpressed[R10LB_0vs0.5_comparison$logFC > 0.6 & R10LB_0vs0.5_comparison$PValue < 0.05] <- "UP"
# if logFC < -0.6 and PValue < 0.05, set as "DOWN"
R10LB_0vs0.5_comparison$diffexpressed[R10LB_0vs0.5_comparison$logFC < -0.6 & R10LB_0vs0.5_comparison$PValue < 0.05] <- "DOWN"

p4 <- ggplot(data=R10LB_0vs0.5_comparison, aes(x=logFC, y=-log10(PValue), col=diffexpressed)) + geom_point() + theme_minimal()

p5 <- p4 + geom_vline(xintercept=c(-0.6, 0.6), col="red") + geom_hline(yintercept=-log10(0.05), col="red")

p6 <- p5 + scale_color_manual(values=c("blue", "black", "red"))
p6 <- p5 + scale_colour_manual(values = mycolors)

R10LB_0vs0.5_comparison$delabel <- NA
R10LB_0vs0.5_comparison$delabel[R10LB_0vs0.5_comparison$diffexpressed != "NO"] <- R10LB_0vs0.5_comparison$gene_symbol[R10LB_0vs0.5_comparison$diffexpressed != "NO"]

ggplot(data=R10LB_0vs0.5_comparison, aes(x=logFC, y=-log10(PValue), col=diffexpressed, label=delabel)) + geom_point() + theme_minimal() +geom_text()

library(ggrepel)
ggplot(data=R10LB_0vs0.5_comparison, aes(x=logFC, y=-log10(PValue), col=diffexpressed, label=delabel)) +geom_point() + theme_minimal() +geom_text_repel() + scale_color_manual(values=c("blue", "black", "red")) + geom_vline(xintercept=c(-0.6, 0.6), col="red") +geom_hline(yintercept=-log10(0.05), col="red")

```

This is the volcano plot I currently have with all of my data:

Can you explain how some of these functions filtered my data (i.e. filterByExpr, estimateDisp, makeContrasts, glmFit, glmLRT)? Why are some of my data points that change from 0 TPM in one condition to some value in another condition showing up on my plot while others are not? Or should I accept that my data analysis is valid or "correct"?

# Answer
edgeR functions do not do any filtering, i.e., they do not remove any genes from the results or from the plots. The only filtering was that which you did explicitly yourself by:

```
keep <- filterByExpr(R10LB_0vs0.5)
R10LB_0vs0.5_filter <- R10LB_0vs0.5[keep, keep.lib.sizes = FALSE]

```

For your data, `filterByExpr` will keep any gene that is expressed at a minimal level in at least 7 samples. You could relax that requirement, to allow genes to be expressed in only 5 or 6 samples say, by reducing the `large.n` argument to `filterByExpr`.
> 1 votes

# Answer
*\[disclaimer: I understand DESeq2, not EdgeR\]*

Most likely, low expression values (as in raw counts, prior to any normalisation such as TPM). Expression that has high dispersion across all samples is less likely to be biologically significant.

It's easier to work out what's going on with an MA plot, with average gene expression on the X axis, and log2 fold change on the Y axis.

I'm used to DESeq2 results, so I'd guess the construction would be something like this (assuming that the 'baseMean' field exists in EdgeR results):

```
R10LB_0vs0.5lrt %>%
    ggplot() +
    aes(y = logFC, col = padj, x=log(baseMean)) +
    geom_point()

```

It's possible that EdgeR doesn't have anything like baseMean in its default results. As a proxy for this, you could use the count totals from the input matrix (e.g. `R10LB_0vs0.5 %>% rowSums() %>% log2()`).
> 0 votes

---
Tags: r, rna-seq, bioconductor, edger, limma
---


# Question
Title: Complete List of IUPAC Chemical Naming Conventions for Formal and Trivial Names?
Two-parted question:

1. Where is a cheatsheet or *complete* list of the "formal" IUPAC naming conventions?
2. Where is a cheatsheet or complete list of the "trivial" or "informal" chemical names?

Does the Wikipedia page have *all* of the naming conventions, or are some missing? Then this page titled "Trivial and semisystematic names retained for naming organic compounds" seems to link to some of the trivial names, but it is in a terrible format for a computer to read (I would like a spreadsheet or some sort of structured data format so I can use these in code) and I'm not sure it's complete. I have searched around on the web but don't seem to find anything that is complete and ideally as structured data. Does anything like this exist? If not structured data, what is the easiest document to read that is complete that I could manually convert into structured data? Also, any of this should ideally be a free resource.

# Answer
> 2 votes
## IUPAC Reports

The International Union of Pure and Applied Chemistry (IUPAC) produces recommendations on nomenclature for chemical compounds. They produce complete reports, and brief guides. However, creating a systematic name from an arbitrary structure is not simple.

It probably depends what 'space' of structures you want to deal with. For simple organic compounds (branched, fused rings, non-bridged) it might be feasible, but there are many edge cases and tricky topics, such as stereochemistry, ligand priority ordering, multiple bridges, and so on.

For example, this is the 'Blue Book' : https://iupac.qmul.ac.uk/BlueBook/ which is a set of PDF docs. Converting these to formal rules would be time consuming, to say the least!

## Informal Names

A resource like ChemSpider provides multiple informal names for compounds, with an API to retrieve data. Seems like (from a quick test) you can get the 'CommonName' field from the details. For example, for a random recordId:

https://api.rsc.org/compounds/v1/records/325/details

you get this return value:

{"id": 325,"commonName": "2-Methoxy-4-vinylphenol"}

so it gives you one common name. Seemingly also for benzene (id 236).

---
Tags: database, cheminformatics
---


# Question
Title: Best methods for clustering heatmap based on regions where the genes are most highly expressed?
I am analysing RNA-Seq data for a set of genes from different human brain regions using data from the Allen Brain Atlas. I have made a heatmap using **GraphPad Prism**. In this heatmap, I show the average RPKM value for each gene of interest for each brain region (the average is taken as there are multiple samples for each brain region).

I am new to creating heatmaps but I know that you can use different methods to cluster the heatmap so that the data is more easy to interpret. I want to cluster the heatmap so that the brain regions which have the highest level of gene expression overall are grouped together.

I have seen online that `R` can also be used to create a heatmap, and I have had previous experience with `R`. Does anyone have any recommendations on any good tutorials (that are not too difficult) on creating heat-maps in `R` as well as clustering the data?

# Answer
> 3 votes
For visualization of expression values it usually makes sense to standardize the data first. The problem with plain expression values is that they differ a lot between genes and as such highly-expressed genes dominate the color scale. As you see in your heatmap, the first and last rows barely show any difference because the color scale is dominated by a single gene CHL1. Therefore the clustering would be due to expression level and not due to differences between columns. For this we standardize the data. Assuming your expression matrix had RPKMs (or similar values) we first log2-transform it and then standardize:

```

mat_log2 <- log2(mat+1)

#/ standardize (Z = Z-scored)
mat_Z <- t(scale(t(mat_log2)))

#/ cluster rows:
hclust_row <- hclust(dist(mat_Z), method="ward.D2")

#/ cluster columns (optional)
hclust_col <- hclust(dist(t(mat_Z)), method="ward.D2")

#/ plot heatmap:
library(ComplexHeatmap)

hm <-
  Heatmap(matrix=mat_Z,
          cluster_rows=hclust_row,
          cluster_columns=hclust_col,
          name="relative expression")
draw(hm)

```

See the ComplexHeatmap documentation for details on the package itself, it is very powerful but also monolithic. This is just an example, any other heatmap package will do fine as well.

Important is to standardize data to put emphasis on differences between columns rather than on absolute expression levels.

---
Tags: r, heatmap
---


# Question
Title: Tutorials on phasing and imputing low-coverage sequencing data
I am new to low-pass whole genome sequencing and have the basic idea of phasing and imputation. I have .vcf file after calling haplotypecaller tools from GATK. After searching the phasing and imputation tools I've encountered some tools like Beagle, Minimac, Shapeit, Glimpse, Eagle. The problem is that I didn't find enough hands-on material on how to process a .vcf from unphased to phased and then to an imputed one.

Can anyone suggests any hands-on tutorial regarding phasing and imputation?

# Answer
This is all relevant for data 0.5x-1x coverage.

Assuming you have genotype likelihoods data, if you want to phase low-coverage data, the most suitable options are GLIMPSE, Beagle4 and QUILT. Of the three, GLIMSPE and QUILT are the most recent and similar in performance. shapeit, eagle and minimac aren't suitable for low-coverage data.

I have found GLIMPSE to be fast and memory efficient, but it requires several different steps; breaking the genome into chunks, imputing, reforming the chunks and then sampling from the haplotypes to phase the data. It has a fairly good tutorial here.

QUILT has fewer stages and faster, but is unlikely to make much of a difference unless you are working on very large datasets. It has a tutorial here.

Which one you use probably doesn't matter a whole lot, but I would probably suggest using QUILT as it has fewer stages and is slightly more recent.
> 3 votes

---
Tags: vcf, imputation, phasing
---


# Question
Title: Pick matching entry from snakemake config table
I'm trying to run a command using snakemake and a TSV config file. Here is my sample config file:

```
sample_id       FQ1_id  FQ2_id
MyCleanName-1       MyFqName-1_R1.fastq.gz        MyFqName-1_R2.fastq.gz
MyCleanName-2       MyFqName-2_R1.fastq.gz        MyFqName-2_R2.fastq.gz

```

Here is my snakefile:

```
import pandas as pd
samples = pd.read_table("samples.tsv").set_index("sample_id", drop = False)

rule all:
    input:
        expand("{sample}_{num}.fastq.gz", sample = samples.sample_id, num = [1,2])

rule my_job:
    input:
        fq1 = expand("{fq1_id}", fq1_id = samples.FQ1_id),
        fq2 = expand("{fq2_id}", fq2_id = samples.FQ2_id),
        logdir = expand("{sample}.logs", sample = samples.sample_id)
    output:
        "{sample}_1.fastq.gz",
        "{sample}_2.fastq.gz",
    message:
        "Running MyCmd on {wildcards.sample}"
    shell:
        """
        MyCmd --sample-name {wildcards.sample} -i1 {input.fq1} -i2 {input.fq2}
        """

rule create_logdir:
    output:
        directory("{sample}.logs")
    message:
        "Creating log directories"
    shell:
        "mkdir -pv {wildcards.sample}.logs/{{MyCmd,MyOtherCmd}}"

```

When I run `snakemake -np MyCleanName-1.fastq.gz`, I get the following:

```
Building DAG of jobs...
Job counts:
        count   jobs
        2       create_logdir
        1       my_job
        3

[Fri Dec 24 16:38:54 2021]
Job 1: Creating log directories

mkdir -pv MyCleanName-1.logs/{MyCmd,MyOtherCmd}

[Fri Dec 24 16:38:54 2021]
Job 2: Creating log directories

mkdir -pv MyCleanName-2.logs/{MyCmd,MyOtherCmd}

[Fri Dec 24 16:38:54 2021]
Job 0: Running MyCmd on MDA-IBC-3

        MyCmd         --sample-name MyCleanName-1 -i1 MyFqName-1_R1.fastq.gz MyFqName-2_R1.fastq.gz -i2 MyFqName-1_R2.fastq.gz MyFqName-2_R2.fastq.gz

Job counts:
        count   jobs
        2       create_logdir
        1       my_job
        3
This was a dry-run (flag -n). The order of jobs does not reflect the order of execution.

```

How can I restrict this so the run only creates log files for the Clean Name being used and uses only the corresponding FASTQs? I feel like I'm missing something trivial here. Should I not be using the config file variable `samples` directly in the rule `my_job`, but somehow wiring up the corresponding `sample_id` and `FQ_id`s on the `all` rule?

# Answer
> 1 votes
Using pointers provided by @bli, I made a few changes and the Snakefile started working as expected.

Instead of using `expand` for the rule's input, I switched to a function that accepts the wildcard sample as input and returns a dict with the actual inputs.

```
 # Could be a one-liner but I wrote it in a verbose
 # manner so it's easier to read
 def get_input_files(in_sample_id):
     fq1 = samples.loc[str(in_sample_id), 'FQ1_id']
     fq2 = samples.loc[str(in_sample_id), 'FQ2_id']
     logdir = "{in_sample}.logs".format(in_sample=str(in_sample_id))
     return_dict = { 'fq1' : fq1, 'fq2' : fq2, 'logdir' : logdir }
     return(return_dict)

```

And the rule `my_job` is now changed (just the `input` part):

```
rule my_job:
    input:
        unpack(get_input_files)
    output:
        "{sample}_1.fastq.gz",
        "{sample}_2.fastq.gz",
...
...

```

Hope this helps someone looking for a similar solution in the future!

---
Tags: snakemake
---


# Question
Title: makeblastdb creating multiple files of unexpectedly large sizes
I have a set of 100 amino acid sequences and I want to perform a BLASTP sesrch against the `refseq_protein` database. Accordingly I had set up the standalone version of BLAST (Version 2.11.0+) and downloaded the refseq\_protein database from NCBI using the following code

```
wget ftp://ftp.ncbi.nlm.nih.gov/refseq/release/complete/*.faa.gz

```

The database gets downloaded as 3027 zipped files containing FASTA sequences. I unzipped all these files and concatenated them into a single file `refseq_protein.faa` (which is around 95 GB in size). Now when I run the following Python code

```
from Bio.Blast.Applications import NcbimakeblastdbCommandline
from Bio.Blast.Applications import NcbiblastpCommandline

cline = NcbimakeblastdbCommandline(dbtype = "prot", input_file = "D:\\refseq_protein.faa", out
= "refseq_protein")

blastp_cline = NcbiblastpCommandline(query = "D:\\DEP_sequences.fasta', db =
"refseq_protein", evalue = 0.01, outfmt = "7 sseqid evalue qcovs pident")

cline()

response = blastp_cline()

```

the `NcbimakeblastdbCommandline` function keeps creating multiple `.phr`, `.pin`, `.psq` etc files which take up a lot of space (In a demo run it had created ~30GB of these files and was still running). I'm afraid this will exhaust the entire space available on my internal hard drive. I'm wondering if there is a way to estimate the total size of the files which `NcbimakeblastdbCommandline` would create. This will help me in deciding whether or not to switch to an external storage to perform the BLASTP search.

I am aware of the fact that pre-formatted refseq\_protein database exists but I'm not sure what value is to be passed in the `db` parameter of the `NcbiblastpCommandline` function, because it asks for the name of the database against which the BLASTP search is to be performed. In the approach that I chose, I had the liberty to set the name of the database.

Any suggestions on how to solve this issue would be appreciated.

# Answer
Welcome to the Bioinformatics site! Hopefully you find the below helpful.

Regarding the creation of a RefSeq protein database locally...as you mentioned, there are pre-made RefSeq BLAST databases available, and looking at the website, the current DB (as of 12/27/2021) is roughly 62.9GB compressed - and is probably ~2x larger decompressed, so I’d guess that the database you are constructing manually will be at least 125GB in size - and likely more.

The above is regarding the technical issue you are currently facing - but I'd back up and look closer at what you are trying to do - align the 100 amino acid sequences to RefSeq. Here are some alternatives that may be -much- easier.

1. BLAST your sequences on the NCBI website. This may not be as satisfying as building out your local tools and databases, but may get the job done more efficiently. This of course depends on whether this is a one-off thing, or something you plan to do frequently, in which case a local database may make sense.
2. Download the premade RefSeq BLAST database. You mentioned you were having issues with how to refer to this database (which may be a useful question in and of itself), but generally an alias file exists within the premade database (entitled `refseq_protein.pal`), and calling this file will allow you to BLAST all the separate chunks. So, the command line might be something like this:

```
blastp -db refseq_protein -query amino_acid_100.fa

```

There are a ton more options you may want to use (try `blastp -help`) , but the above should work.

3. Use the `-remote` option in blastp to run from the command line against the NCBI server. The below line should perform the BLAST alignment against the NCBI refseq protein database.

```
blastp -remote -db refseq_protein -query amino_acid_100.fa

```
> 2 votes

---
Tags: fasta, python, refseq, blastp, makeblastdb
---


# Question
Title: Pros and cons between genomes and transcriptomes?
What are the major differences between the use cases of a shot-gun genome and an mRNA capture transcriptome? Especially when it comes to downstream analyses such as looking for orthology, and selection analysis.

From what I can visualize in my head, I can come up with the following points:

1. The genome has the structural data, information on chromosomes, and haplotypes while the transcriptome lacks this information.
2. The genome records the SNPs more directly than transcriptomes. I might be thinking of this wrong but the genomes would have the actual bases and we can then compare the bases between the genomes but for transcriptomes differences in expression would show up as presence/absence variations that have no way of being traced back?
3. The transcriptomes should be easier to work with when it comes to computational requirements since they tend to be 25 to 50 times smaller than genomes.

I could be completely wrong about the points I listed and could be missing something major. What do you guys think?

# Answer
It really depends on your use case. If this is human data, where we already have a fairly polished reference genome, then genome sequencing adds value for the genetic variation, e.g. rare disease genetics, GWAS or cancer genomics. Transcriptome sequencing allow exploration of tissue (or cell type) specific patterns and RNA biomarkers etc. Calling variants from RNA can be interesting in the context of splicing and RNA editing.

For a non-model organism where there isn't currently a reference genome, then the genome creates a reference for future studies to map genomic/transcriptomic reads to, so it's usual to do the genome first.

Since you mention ortholgy and selection studies, well traditional comparative genomics studies have done these analyses from the genome only, but the transcriptome can add value depending on your questions.

Finally, I wouldn't consider computational requirements as you stated: it depends on sequencing depth/coverage and other factors, and is likely not to drive the choice between the two assays
> 0 votes

---
Tags: transcriptome, genome-sequencing, orthologues
---


# Question
Title: Help with MinION sequencing data species identification
Hi I'm new to bioinformatics and have just completed my first run on the MinION (long read sequencing Oxford Nanopore Technologies). I was hoping someone could direct me towards R packages, workflow, tutorials or guides that will help me identify species that are present in my sample mainly for fungi as my amplicons had ITS1 and ITS2 regions.

I have been watching Youtube tutorials on RStudio for over a week and I feel like I'm no closer than I was when I started. I have both fastq and fast5 files and I tried to run the fastq files using EPI2ME, however, the database is lacking in fungi. I have tried using NCBI's BLAST but I can only run \<5kb reads and most of my sequencing reads are \>6kb. Also, I've tried using Galaxy to convert my fastq files to fasta but I'm not sure how to run my long reads for species identification using Galaxy or if that is even possible.

Any help would be greatly appreciated. Thanks.

# Answer
There are a number of useful metagenomics tools out there, including Centrifuge and Kraken2, as @gringer has mentioned.

Also, FWIW, the WIMP tool from ONT, a Centrifuge based metageonomic classifier, does include fungus sequences derived from RefSeq, though it may not include the species you are looking for.

If you want to use NCBI BLAST, you could split your long reads into pieces that are \<5kb (there are a number of tools out there for FASTQ manipulation - my current favorite tool for this is seqkit, where you could do something like this):

`seqkit sliding -s 4900 -W 4900 test1.fastq`

The above will split your FASTQ reads into chunks of 4.9kb. You can adjust the `-s` parameter if you instead want overlapping chunks.

From a GUI standpoint, you can use Galaxy for the above by making a workflow with the following steps:

1. Trim Sequences (enter reasonable numbers here, to trim each of your reads to \<5kb)
2. FASTQ to FASTA converter
3. NCBI BLAST+

Another possibility is QIIME2, which can be used with the UNITE fungal training set. See recent reference (looking at the ITS1/ITS2 regions), UNITE, and QIIME2.

Long term, you may want to look into using the command line, as it is much more powerful and customizable, but if BLAST gives you reasonable results for your reads, the above process in Galaxy may work for you.
> 1 votes

# Answer
I've found Centrifuge and Kraken2 to be useful for species classification. There are pre-prepared databases that include Fungi, e.g. see here for the Kraken 2 "PF" database, which is the Standard database plus protozoa & fungi.

Can't really provide much more help than that with more information about what you've done.

You'll need a computer with around 64 GB of memory in order to be able to run these programs, and it's probably not the best idea to jump right into command line compilations when starting out... but it looks like you're up for a challenge given that your first MinION run has involved jumping into the deep end with a fungi metagenomic project.
> 0 votes

---
Tags: sequence-alignment, sequence-analysis, minion, galaxy
---


# Question
Title: Nextflow HISAT2 command exit status 255 with prebuild index
I have a problem with my nextflow pipeline.

My workflow looks like this:

1. Uncompress the hisat2 genotpye index grch38\_genome.tar.gz from https://genome-idx.s3.amazonaws.com/hisat/grch38\_genome.tar.gz
2. Align .fastq.gz files from /https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR1806626 and https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR1806627

However, in process HISAT2 I get command exit status 255:

> Command error: "grch38" does not exist

The command executed is:

> hisat2 --very-sensitive --no-spliced-alignment -x grch38 -U SRR1806626.fastq.gz \> SRR1806626.fastq.sam

Can someone please explain to me what I'm doing wrong?

My script.nf file for reproducing the error looks like this:

```
#!/usr/bin/env nextflow

nextflow.enable.dsl=2

params.reads = "$baseDir/data/*{1,2}*.fastq.gz" //https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR1806626, https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR1806627
params.hisat2_index = "$baseDir/data/*_genome.tar.gz" // https://genome-idx.s3.amazonaws.com/hisat/grch38_genome.tar.gz

log.info """\
         reads:        ${params.reads}
         hisat2_index: ${params.hisat2_index}
         """
         .stripIndent()

process UNCOMPRESS_GENOTYPE_INDEX {

  tag "$hisat2_index_ch.baseName"

  input:
  path(hisat2_index_ch)

  output:
  path('*')

  script:
  """
  tar xvzf $hisat2_index_ch
  """

}

process HISAT2 {

  tag "$reads.baseName"

  input:
  path(hisat2_indices)
  path(reads)

  output:
  path('*.sam')

  script:
  """
  hisat2 --very-sensitive \
         --no-spliced-alignment \
         -x $hisat2_indices \
         -U $reads \
         > ${reads.baseName}.sam
  """

}

workflow {

  Channel
    .fromPath( params.reads )
    .ifEmpty { exit 1, "reads was empty - no input files supplied" }
    .set { reads_ch }

  Channel.fromPath( params.hisat2_index )
    .ifEmpty { exit 1, "hisat2_index was empty - no input file supplied" }
    .set { hisat2_index_ch }

  UNCOMPRESS_GENOTYPE_INDEX(hisat2_index_ch)

  HISAT2(UNCOMPRESS_GENOTYPE_INDEX.out, reads_ch)

  //Command error:
  //(ERR): "grch38" does not exist
  //Exiting now ...

}

```

Kind regards, Thomas

# Answer
> 0 votes
The problem is probably that the tar process creates a folder with the genome files but the input for `-x` is not the folder with the index files but the basename of the index files itself.

---
Tags: rna-seq, nextflow, rna-alignment, errors
---


title: Use of "baseMean" in DeSEQ2

text:
# Question
Title: Use of "baseMean" in DeSEQ2
The basemean is described as the "mean of normalized counts of all samples". My question is, how would we interpret differences in basemean between genes, pertaining to reliability and such? For example, if two genes that have similar p-values and l2fchanges, but one has a much higher basemean, would that gene's data be more dependable?

Sorry if this is a trivial question!

# Answer
I could write a long story about this but actually I will just link this excellent resource from HBCtraining which dissects the individual steps DESeq2 does, starting from normalization over dispersion estimation (this is where the baseMean comes into play) over model fitting and testing. I think it will clarify the role of the baseMean in DESeq2.
> 2 votes

---
Tags: deseq2
---


# Question
Title: Interpreting imputation result from GLIMPSE
I'm following this tutorial of GLIMPSE for learning. I was expecting some extra SNPS coming from the 1000 genome reference in the resulting .vcf file. Though I understand the phasing in the output resulting .vcf, what I don't understand is the imputed result. Apparently the .vcf before and after has the same number of SNPs. I ran the scripts from here. I'm including the vcf files before and after imputation.

Can anyone explain me where the imputed SNPs are in the output v.cf file. If my concept is wrong can anyone tell me what should we expect in result.

# Answer
It seems that you didn't use the following flag:

```
 --impute-reference-only-variants   Allows imputation at variants only present
                                     in the reference panel. The use of this
                                     option is intended only to allow
                                     imputation at sporadic missing variants.
                                     If the number of missing variants is
                                     non-sporadic, please re-run the genotype
                                     likelihood computation at all reference
                                     variants and avoid using this option,
                                     since data from the reads should be used.
                                     A warning is thrown if reference-only
                                     variants are found.

```

You were expecting additional SNPs which were present in the reference panel to be imputed into your targets, but you didn't include the flag `--impute-reference-only-variants`, so only individual missing genotypes were imputed. Individual missing genotypes `./.` are always imputed if they are present in both the target and reference, but if the SNP is present in the reference and not the target, it won't be unless you use that flag.
> 2 votes

---
Tags: vcf, imputation, phasing, 1000genomes
---


# Question
Title: RNA-Seq Data Heatmap: Is it necessary to do a log2 transformation of RPKM values before doing the Z-score standardisation?
I am making a heatmap using RNA-Seq data in R. The heatmap shows gene expression values (RPKM) in different brain regions. I have the following code:

```
library(tidyverse)

library(pheatmap)

library(matrixStats)

read_csv("prenatal_heatmap_data.csv") -> all_data

all_data %>%

   column_to_rownames("Brain Region") -> heatmap_data

heatmap_data %>%

   pheatmap()

```

Which generates the following heatmap:

I want to do a log2 transformation of the RPKM values so that I can do the Z-score standardisation of the data. I have log2-transformed the values using the following code:

```
heatmap_data %>%

   log2() -> heatmap_data_log2

heatmap_data_log2 %>%

   pheatmap()

```

However, when I try to create a heatmap using the log2-transformed values I get the following error:

I have looked at the log2-transformed data and I know that the reason I am getting this error is because in the original dataset, some of the RPKM values were 0. When these 0 values get log-transformed they become -Inf:

^ Original Data

^ log2 transformed data

I am not sure how to overcome this issue. I was wondering if it is necessary to do the log2 transformation of the RPKM values before doing the Z-score standardisation? I have seen that it is conventional to do the log2-transformation of the data.

Any advice is appreciated.

# Answer
> 1 votes
Add a pseudocount of one as in the code suggestion I gave you earlier: Best methods for clustering heatmap based on regions where the genes are most highly expressed?

Log of zero is not defined.

---
Tags: r, heatmap
---