# Question
Title: Which sequence alignment tools support codon alignment?
Sometimes it useful to perform a nucleotide protein coding gene sequence alignment based on codons, not on individual nucleotides. For example for further codon model analysis it is important to have full codons.

A widely used approach here is to perform a protein sequence alignment first and then impose this alignment to the nucleotide sequences using PAL2NAL, CodonAlign or something similar.

This is how transAlign or GUIDANCE (in codon mode) work.

The problem here is that you are discarding part of the information which could be potentially used for the sequence alignment. E.g. if you have slowly evolving low-complexity region adjacent to a quickly evolving one, the amino acid induced alignment could be wrong, while incorporating nucleotide sequence potentially allows to make the alignment more accurate.

I'm aware of two programs which can do true codon alignment. First, PRANK has a dedicated codon model, but it is rather slow and using it is overkill for certain problems. Second, Sequence Manipulation Suite can perform codon alignments, but only for a pair of sequences; also it's javascript based, therefore it is hard to run it for a large number of sequences.

Can you recommend any software for multiple codon sequence alignment? Preferably available for offline use.

# Answer
I don't know of any transcript-to-transcript aligners that are able to do this, but LAST can align transcript queries to protein reference sequences using a specified frameshift cost. Here's the specific documentation for that option:

> -F COST
>
> Align DNA queries to protein reference sequences, using the specified frameshift cost. A value of 15 seems to be reasonable. (As a special case, -F0 means DNA-versus-protein alignment without frameshifts, which is faster.) The output looks like this:
>
> ```
> a score=108 s prot 2  40 + 649
> FLLQAVKLQDP-STPHQIVPSP-VSDLIATHTLCPRMKYQDD s dna  8 117 + 999
> FFLQ-IKLWDP\STPH*IVSSP/PSDLISAHTLCPRMKSQDN
>
> ```
>
> The \ indicates a forward shift by one nucleotide, and the / indicates a reverse shift by one nucleotide. The * indicates a stop codon. The same alignment in tabular format looks like this:
>
> ```
> 108 prot 2 40 + 649 dna 8 117 + 999 4,1:0,6,0:1,10,0:-1,19
>
> ```
>
> The "-1" indicates the reverse frameshift.

I sent an email to the LAST mailing list about adding a frameshift penalty for transcript-to-transcript matching; I've been pleasantly surprised with the requested features that Martin Frith has added to LAST in the past. Unfortunately, in this case the problem is too difficult to sort out due to all the possible combinations that could happen, so it's unlikely to be implemented in LAST in the forseeable future (unless someone else writes that code).
> 4 votes

# Answer
Try MACSE v2 (https://academic.oup.com/mbe/article/35/10/2582/5079334) will align multiple protein-coding nucleotide sequences based on their amino acid translation while allowing for the occurrence of frameshifts
> 3 votes

# Answer
Virulign does this for virus sequences, the publication is available here and github address here
> 2 votes

---
Tags: sequence-alignment, covid-19, codon, sars-cov-2, virulign
---


# Question
Title: Output the linear predictor from a stratified cox model?
I have fitted a cox model on a pooled dataset of multiple studies, say studies A, B, and C. As these studies have different baseline hazard functions, I stratified for 'study' in the model.

Like this:

```
df$SurvObj <- with(df, Surv(event_rd, event == 1))

fit <- coxph(SurvObj ~ cov1 + cov2 + cov3 + cov4 + strata(study),data=df)

```

Now, I want to assess whether I can use the beta coefficients from the above model to predict the event probability in a new study, study D. The baseline hazard in study D is different from those in studies A, B, and C.

When I have

```
predict(fit,type="lp",newdata=studyD,reference="strata")

```

I get this error:

> Error in model.frame.default(data = studyD, formula = ~hba1c + sbp + : factor strata(study) has new levels study=4

Why does R require that I match the strata of study D with those in study A, B, and/or C when I only want to output the linear predictor (`type="lp"`)? As the linear predictor is independent of the baseline hazard function in this case.

I suppose it should be possible to extract the linear predictor for individuals in study D and then manually calculate the event probability using the baseline hazard of study D.

Does anybody know how to do this?
I could not find any question like this so maybe I am thinking the wrong way.

# Answer
Unlike wikipedia, the linear predictor (LP) in `coxph` is made dependent from the baseline hazard function $\lambda$. For example, if a model has two predictors `x1`, `x2` stratified by `sex`, the linear predictor would be like: $$LP = log ( \lambda \_{sex}) + \beta \_{1}x \_{1}+\beta \_{2}x \_{2}$$

Where the $log(\lambda \_{sex})$ can be obtained by `predict()` with `x1=x2=0`.

We could calculate the LP without strata information (sex) by hand using the $\beta \_{1} \\, \beta \_{2}$ provided by `coef(fit)` and a customised $\lambda \_{unkown \\, sex}$.

Below is a simulation comparing the results from `predict()` and by hand `coef(fit)`.

```
## Toy data
df <- list(time=c(4,3,1,1,2,2,3),
          status=c(1,1,1,0,1,1,0),
          x1=c(0,2,1,1,1,0,0),
          x2=1:7,
          sex=c(0,0,0,0,1,1,1))

## Fit a stratified model
fit <- coxph(Surv(time, status) ~ x1 + x2 + strata(sex), df)
coef(fit)
#         x1          x2
# 0.79451881 -0.01917633

## Baseline linear predictor
lambda_sex1 <- predict(fit, newdata=list(x1=0, x2=0, sex=0))
# -0.746578
lambda_sex2 <- predict(fit, newdata=list(x1=0, x2=0, sex=1))
# -0.1497816

## by predict function
predict(fit, newdata=list(x1=1, x2=2, sex=0))
# 0.009588167
predict(fit, newdata=list(x1=1, x2=2, sex=1))
# 0.6063845

## by hand
lambda_sex1 + coef(fit)[1]*1 + coef(fit)[2]*2
# 0.009588167
lambda_sex2 + coef(fit)[1]*1 + coef(fit)[2]*2
# 0.6063845

```

Finally, to calculate "a" linear predictor for an unknown `sex`:

```
# prevented by predict function
predict(fit, newdata=list(x1=1, x2=2, sex=2))
#Error in model.frame.default(data = list(x1 = 1, x2 = 2, sex = 2), formula = ~x1 +  :
# factor strata(sex) has new level sex=2

# by hand, assuming lambda for sex 2 is 0
0 + coef(fit)[1]*1 + coef(fit)[2]*2
# 0.7561661

```
> 3 votes

---
Tags: r
---


# Question
Title: Download multiple fastq files using fastq-dump
I want to download the following fastq files at the same time in Salmon:

```
 - SRR10611214
 - SRR10611215
 - SRR10611215
 - SRR10611216
 - SRR10611217

```

Is there a way to do this using a bash for loop or fastq-dump? or prefetch

# Answer
> 4 votes
A sample code is given in the salmon documentation as follows. Source

```
#!/bin/bash
mkdir data
cd data
for i in `seq 25 40`;
do
  mkdir DRR0161${i};
  cd DRR0161${i};
  wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/DRR016/DRR0161${i}/DRR0161${i}_1.fastq.gz;
  wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/DRR016/DRR0161${i}/DRR0161${i}_2.fastq.gz;
  cd ..;
done
cd ..

```

This could be modified as follows.

```
#!/bin/bash
mkdir data
cd data
for i in `seq 14 17`;
do
  wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR106/0${i}/SRR106112${i}/SRR106112${i}_1.fastq.gz;
  wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR106/0${i}/SRR106112${i}/SRR106112${i}_2.fastq.gz;
done
cd ..

```

You can save the code as a shell script and run it from the linux terminal. For example `bash download_fastq.sh`

# Answer
> 1 votes
Make a list.txt file containing a single column of SRA numbers to download.

then:

`for i in $(cat list.txt); do echo $i; date; fasterq-dump -S $i; done`

It works well to use NCBI's web interface to find SRA samples of interest, download and open findings in Excel, then copy single column containing SRA numbers and paste into list.txt using document editor such as vim.

After downloading including the "R" can be nice:

`for i in *_1.fastq; do mv $i ${i%_1.fastq}_R1.fastq; done`

`for i in *_2.fastq; do mv $i ${i%_2.fastq}_R1.fastq; done`

and zip:

`pigz *fastq`

If needed a conda option for downloading fasterq-dump:

`conda install -c bioconda sra-tools`

# Answer
> 0 votes
You can use `parallel`.

```
parallel -j 3 fastq-dump {} ::: SRR10611214 SRR10611215 SRR10611215 SRR10611216 SRR10611217

```

The option `-j` says how many jobs should maximal run parallel. So in this case maximal 3 identifier would be handled at the same time.

How many jobs you can run parallel depends on your machine.

You can also take a look at parallel-fastq-dump.

---
Tags: rna-seq, fastq, data-download, salmon
---


# Question
Title: How to fit all genes (labels) in chromosome ideogram plot made by RCircos package?
I am using Rcircos to make a chromosome ideogram plot for my gene list (n = 45). However, I am getting this error:

```
 Not all labels will be plotted.
 Type RCircos.Get.Gene.Name.Plot.Parameters()
 to see the number of labels for each chromosome.

```

I am using the below-given script and what to know is there any way to fit all my genes in the chromosome ideogram plot

```
 library(RCircos)
 data(UCSC.HG19.Human.CytoBandIdeogram)
 RCircos.Set.Core.Components(
 cyto.info=UCSC.HG19.Human.CytoBandIdeogram,
 chr.exclude = NULL,
 tracks.inside = 10, tracks.outside = 0)

 #Making a plot with RCircos
 out.file <- "RCircos_IF.pdf";
 pdf(file=out.file, height=8, width=8, compress=TRUE);
 RCircos.Set.Plot.Area();
 par(mai=c(0.25, 0.25, 0.25, 0.25));
 plot.new();
 plot.window(c(-2.5,2.5), c(-2.5, 2.5));

 # plot chromosome ideogram
 RCircos.Set.Plot.Area();
 RCircos.Chromosome.Ideogram.Plot()
 gene.lable.data = read.table ("IF.txt", sep="\t", header=T, row.names = 1) # IF.txt is my input file
 RCircos.Gene.Connector.Plot(gene.lable.data, track.num = 1, side = "in");
 track.num <- 2;
 RCircos.Gene.Name.Plot(gene.lable.data,name.col = 4, track.num)

 dev.off()

```

My input file looks like

```
    Chromosome  chromStart  chromEnd    Gene.name
    1   chr1    150363091   150476566   RPRD2
    2   chr1    150549369   150560937   ADAMTSL4
    3   chr1    91949371    92014426    BRDT
    4   chr1    31365625    31376850    FABP3
    5   chr1    150960583   150975004   CERS2

```

# Answer
> 2 votes
To change the maximum number of genes in `RCircos`, we could modified the `char.width` according to this link. The actual value in a RCircos session could be modified with get and reset methods for plot parameters.

```
params <- RCircos.Get.Plot.Parameters() #$char.width
params$char.width <- 100 #default 500
RCircos.Reset.Plot.Parameters(params)

#the maxLabels are updated accordingly
RCircos.Get.Gene.Name.Plot.Parameters()

```

# Answer
> 0 votes
Just want to add that apart from params$char.width, we can also use

```
params$text.size
params$text.size <- 0.2 #default 0.4
params$base.per.unit < 50000 #imp #default value = 30000

```

---
Tags: r, circos
---


# Question
Title: can I download DESeq2 in R 3.6.3 in Linux MInt?
I am trying to download DESeq2 in R 3.6.3. Is it possible? This is printed when I am trying: Warning in install.packages : package ‘DESeq2’ is not available (for R version 3.6.3) Thank you!

# Answer
> 1 votes
Bioconductor is a package management tool for R that allows you to install the appropriate package versions for your R version.

For R 3.6.3 you'll need Bioconductor 3.10:

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

```

Then you install DESeq2:

```
BiocManager::install("DESeq2")

```

# Answer
> 0 votes
these errors are generated when I am doing what PPK suggested: Bioconductor version 3.10 (BiocManager 1.30.10), R 3.6.3 (2020-02-29) Installation path not writeable, unable to update packages: nlme, spatial, survival

> BiocManager::install("DESeq2") Bioconductor version 3.10 (BiocManager 1.30.10), R 3.6.3 (2020-02-29) Installing package(s) 'DESeq2' Warning: dependency ‘XML’ is not available also installing the dependencies ‘annotate’, ‘genefilter’, ‘geneplotter’

# trying URL 'https://bioconductor.org/packages/3.10/bioc/src/contrib/annotate\_1.64.0.tar.gz' Content type 'application/x-gzip' length 1873132 bytes (1.8 MB)

downloaded 1.8 MB

# trying URL 'https://bioconductor.org/packages/3.10/bioc/src/contrib/genefilter\_1.68.0.tar.gz' Content type 'application/x-gzip' length 1421076 bytes (1.4 MB)

downloaded 1.4 MB

# trying URL 'https://bioconductor.org/packages/3.10/bioc/src/contrib/geneplotter\_1.64.0.tar.gz' Content type 'application/x-gzip' length 1431450 bytes (1.4 MB)

downloaded 1.4 MB

# trying URL 'https://bioconductor.org/packages/3.10/bioc/src/contrib/DESeq2\_1.26.0.tar.gz' Content type 'application/x-gzip' length 2060061 bytes (2.0 MB)

downloaded 2.0 MB

ERROR: dependency ‘XML’ is not available for package ‘annotate’

* removing ‘/home/nikolena/R/x86\_64-pc-linux-gnu-library/3.6/annotate’ ERROR: dependency ‘annotate’ is not available for package ‘genefilter’
* removing ‘/home/nikolena/R/x86\_64-pc-linux-gnu-library/3.6/genefilter’ ERROR: dependency ‘annotate’ is not available for package ‘geneplotter’
* removing ‘/home/nikolena/R/x86\_64-pc-linux-gnu-library/3.6/geneplotter’ ERROR: dependencies ‘genefilter’, ‘geneplotter’ are not available for package ‘DESeq2’
* removing ‘/home/nikolena/R/x86\_64-pc-linux-gnu-library/3.6/DESeq2’

The downloaded source packages are in ‘/tmp/RtmphLXrZb/downloaded\_packages’ Installation path not writeable, unable to update packages: nlme, spatial, survival Warning messages: 1: In install.packages(...) : installation of package ‘annotate’ had non-zero exit status 2: In install.packages(...) : installation of package ‘genefilter’ had non-zero exit status 3: In install.packages(...) : installation of package ‘geneplotter’ had non-zero exit status 4: In install.packages(...) : installation of package ‘DESeq2’ had non-zero exit status

I have tried to instal XML and the other dependencies without success..

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


title: Changing active.ident in Seurat

text:
# Question
Title: Changing active.ident in Seurat
Im trying to change the active.ident to another column in metadata but this error keeps popping up! I recently upgraded to R version 4.0.2 from 3.6.1 The older version was working but the new one isn't. I don't understand how to fix it?

```
M@active.ident <- M@meta.data$status
Error in (function (cl, name, valueClass)  :
assignment of an object of class “character” is not valid for @‘active.ident’ in an object of class
“Seurat”; is(value, "factor") is not TRUE

head(M@meta.data)
        orig.ident nCount_RNA nFeature_RNA percent.mito   status  RNA_snn_res.0.6 seurat_clusters
BC03_03       BC03   998414.7         3356   0.07869098 nonTumor                0               0
BC03_09       BC03   999853.6         7616   0.12964336  primary                0               0
BC03_13       BC03   998206.5         8457   0.17672205  primary                0               0
BC03_17       BC03   999535.1         7494   0.11323008  primary                0               0

```

I want to run the UMAP with the clusters named according to the status and not seurat clusters.However, the code isnt working..

# Answer
`M <- SetIdent(M, value = "status")`

or more explicitly

`M <- SetIdent(M, value = M@meta.data$status)`

You can also use the `group.by` argument of `UMAPPlot()` or other plotting functions from `Seurat` for that matter.
> 2 votes

---
Tags: r, scrnaseq, seurat, sequence-annotation, clustering
---


# Question
Title: calculateQCMetrics defunct, how to calculate the quality metrics by perCellQCMetrics
I am trying to follow a tutorial from Sanger institute (from May 2019) on analysis of single cell RNA Seq data. They use calculateQCMetric function to calculate the quality metrics, but I am getting an error message that calculateQCMetric is a deprecated (isSpike function is also deprecated). I want to calculate the quality metrics. Can anyone provide code for getting QCMetric step? I attached the previous defunct code below:

```
umi <- calculateQCMetrics(
   umi,
   feature_controls = list(
       ERCC = isSpike(umi, "ERCC"),
        MT = isSpike(umi, "MT")
   )

```

And I use StupidWolf's code from this question isSpike function in SingleCellExperiment package is deprecated? to instead isSpike function. I also attached below:

```
Spikein_names = grep("^ERCC-",rownames(molecules),value=TRUE)
SpikeIn = molecules[Spikein_names,]

g = genes(EnsDb.Hsapiens.v86)
MTgene_names = g[seqnames(g)=="MT"]$gene_id
MTgenes = molecules[rownames(molecules) %in% MTgene_names,]
Have a vector of genes to keep:

keep = setdiff(rownames(molecules),c(SpikeIn,MTgenes))
Then make the single cell object:

sce <- SingleCellExperiment(
assays=list(counts=as.matrix(molecules[keep,])),colData=annotation)
Use altExp to slot in the stuff:

altExp(sce, "spike-in") <- SummarizedExperiment(SpikeIn)
altExp(sce, "MTgenes") <- SummarizedExperiment(MTgenes)

```

I try the below code, which did not work.

```
> perCellQCMetrics(sce,
 subsets = list(Spikein_names,MTgene_names),
 flatten = TRUE,
 exprs_values = "counts"，
 use_altexps = TRUE)

```

Thanks in advance.

# Answer
Ok I have figured it out and I am sharing this in case anyone can get help. I used

```
sce <-addPerCellQC(sce)

```

which calculated the QC for all the altexps also and adds them to the colData.
> 1 votes

# Answer
It would be helpful if you can provide the error message. Guessing from your codes, it's possible that the alternative experiments weren't set correctly. My suggestion:

```
altExp(sce, "spike-in") <- SummarizedExperiment(assays=list(counts=as.matrix(SpikeIn)))
altExp(sce, "MTgenes") <- SummarizedExperiment(assays=list(counts=as.matrix(MTgenes)))

```

If you're using altExp, the subsets argument in perCellQCMetrics is unnecessary in this case.

```
perCellQCMetrics(sce,
exprs_values="counts",
use_altexps=TRUE)

```
> 0 votes

---
Tags: r, scrnaseq, bioconductor
---


# Question
Title: mBED has time complexity $O(n \lg n)$, claimed by Clustal Omega paper, why?
In the Clustal Omega paper, it says the following about mBED ($N$ is the number of input sequences, I assume):

> we use a modified version of mBed (Blackshields et al , 2010), which has complexity of $O(N \log N)$

But in the mBED paper it says the following (section 2):

> A number of $t$ sequences are initially sampled from the input dataset $X .$ Following the LLR algorithm, this value is set by default to $t=\left(\log \_{2} N\right)^{2} .$

> After the seed sequences in $R$ have been chosen, all sequences in the input data are associated with a $t$-dimensional vector. This is done simply by computing the distances from all sequences to each of the $t$ seeds.

Then my conclusion is that this would lead to a complexity of $O(N (\log N)^2)$ because for each sequence (in total $N$ of them), distances to $t = (\log N)^2$ sequences must be calculated. Where did I do wrong?

# Answer
Turns out that I was dumb: I misread.

The Clustal Omega paper (in the supplementary materials) chose $t = \log N$, which I think originates from this claim from this part of the Linial et al. paper:

> In random polynomial time $X$ may be embedded in $l\_{2}^{O(\log n)}$ with distortion $\leq(1+\epsilon) \cdot c\_{2}(X)$ (for any $\epsilon\>0)$

Several questions still remain though:

1. The newest implementation of Clustal Omega (that I downloaded) has source code that instead chooses $t = (\log N)^2$, why changed?
2. The original mBED paper chooses $t = (\log N)^2$, but it was not justified on how they chose this number instead of $\log N$
> 1 votes

---
Tags: sequence-alignment
---


# Question
Title: Samtools Index: Chromosome Blocks not Continuous
I am working with short-read whole-genome sequences from the NCBI's SRA. I have aligned and sorted all of my short-read sequences and am attempting to index each sequence into .bai format using samtools index, but am running into a couple of errors.

I unpacked the original .sra files in the following manner:

```
fastq-dump --defline-seq '@$sn[_$rn]/$ri' --split-files SRR6509138.sra

```

I then aligned each fastq paired-end read, removed duplicate reads, and converted each to bam format like so:

```
bwa mem xiphophorus_birchmanni_10x_12Sep2018_yDAA6.fasta SRR6509136_1.fastq SRR6509136_2.fastq | samblaster -e -r| samtools view -Sb - > blasted_SRR6509136.bam

```

Some of my files did not come as paired-end reads. For those files, I used the --ignoreUnmated flag on samblaster.

I then sorted each bam file like so:

```
samtools sort blasted_SRR6649368.bam -o sorted_SRR6649368.bam -n

```

I am attempting to obtain a bai file for each bam using the following command:

```
samtools index sorted_SRR6649368.bam sorted_SRR6649368.bam.bai

```

This is the error I run into for the unpaired reads:

```
[E::hts_idx_push] Chromosome blocks not continuous
[E::sam_index] Read 'HWI-ST387:164:D0CJWACXX:4:1101:1129:13264' with ref_name='ScdB1pO_646;HRSCAF=880', ref_length=33092979, flags=16, pos=23260108 cannot be indexed
samtools index: failed to create index for "sorted_SRR791885.bam": No such file or directory

```

This is the error I run into for the paired reads:

```
[E::hts_idx_push] Unsorted positions on sequence #79: 11159020 followed by 11158717
samtools index: failed to create index for "sorted_SRR6649368.bam"

```

Can anyone help me figure out why this is happening?

# Answer
```
samtools sort blasted_SRR6649368.bam -o sorted_SRR6649368.bam -n

```

These error messages indicate that the reads are not sorted by coordinate — in particular, that the reads mapped to *ScdB1pO\_646;HRSCAF=880* are not all together, and that the reads at positions 11159020 \< 11158717 are not sorted by position on their chromosome.

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

---
Tags: genome, samtools, bwa, indexing
---


# Question
Title: How to install genometools with anaconda?
I tied to install genometools with anaconda:

```
conda install -c bioconda genometools

```

Though, after running that successfully the command `gt` is not available in my conda environment.

> The program 'gt' is currently not installed. You can install it by typing: sudo apt install genometools

While conda tells me the program is installed. Are there multiple packages called genometools? What is the difference and can I install the package providing `gt` with anaconda at all? There is also a package called `genometools-genometools`, but it only has two downloads. So, I am not sure that is what I need.

# Answer
> 2 votes
`genometools-genometools` seem to be the current version. By installing the `genometools` package I received the same error like you. I guess it's because the binary isn't copy to `bin` during installation.

`genometools-genometools` seems to work fine. And it has `2k` (2000) and not `2` :)

# Answer
> 0 votes
genometools-genometools worked for me, while 2 other options (genometools and -c leylabmpi genometools) did not.

---
Tags: genometools
---


# Question
Title: Basic question about finding mRNA sequence in transcriptome
I am a computer scientist just starting out in bioinformatics topics, and would appreciate any guidance that can be given here:

I have an mRNA sequence\- an isoform - whose length is about 4000 base pairs (bp) - the origin is provided in the link.

I also have a bam file/fasta file of the entire transcriptome of a kidney cell (see Data Access tab).

My objective is to quantify mRNA expression levels of the mRNA sequence in the transcriptome.

How does one get started? Should I use the original format (TenX bam, about 10 GB) or SRA archive data (fasta, about 3.5 GB)? Why are the files so different in size?

I have a MacBook with `samtools` installed and 16GB of RAM. Is this enough to process the data and what basic commands should be used to quantify the 4000 bp sequence overlaps in the transcriptome?

# Answer
Basically you want to know the VHL Gene expression in that 10x dataset. They have provided the barcode, matrix and feature files on GEO GSE131685 as Supplementary file. You could use those files and get the gene expression matrix as described on 10x website.

```
library(Matrix)

matrix_dir = "/opt/sample345/outs/filtered_feature_bc_matrix/"
barcode.path <- paste0(matrix_dir, "barcodes.tsv.gz")
features.path <- paste0(matrix_dir, "features.tsv.gz")
matrix.path <- paste0(matrix_dir, "matrix.mtx.gz")

mat <- readMM(file = matrix.path)

feature.names = read.delim(features.path,
                           header = FALSE,
                           stringsAsFactors = FALSE)

barcode.names = read.delim(barcode.path,
                           header = FALSE,
                           stringsAsFactors = FALSE)

colnames(mat) = barcode.names$V1
rownames(mat) = feature.names$V1

```

This mat can be further loaded into Seurat for any downstream analysis. All this can be done on a laptop.

If you want to re-quantify the gene expression for whatever the reason, you need to re-run the cellRanger pipeline using GTF of your interest (it needs some preparation), which may not work on your laptop. It requires large memory and long running times, ideally to submit as job on cluster.
> 2 votes

# Answer
> Should I use the original format (TenX bam, about 10 GB) or SRA archive data (fasta, about 3.5 GB)? Why are the files so different in size?

I strongly recommend you hold your horses; how do you think you can analyze something if you don't even know what you are looking at?

Since this is a "TenX" bam, you can probably do the job with samtools and grep. Looking up 10XGenomics bam tags might be helpful. You might need to look up SAM flags too.

> My objective is to quantify mRNA expression levels of the mRNA sequence in the transcriptome.

You understand that this isn't really a thing that's possible? With spike-ins, maybe, but you don't have that. All you can do is work out the proportions of reads belonging to that transcript (Or gene, really, unless you want to put in a great deal more work, beyond I believe what 10XGenomics does). Which I guess you want to do by lumping all the cells together, and just pretending there is no QC to do?
> 1 votes

---
Tags: sam, fasta, samtools, transcriptome
---


# Question
Title: Split fasta file based on groups in header information and output as separate files
I have a fasta file containing the sequence of a gene across different species. In total there are around 900 samples and 12 species. (Each sequences is over multiple lines and longer than 100bp.)

My fasta file looks like:

```
>Species-1-samplenameA
CTATCCTTAAACGCATATCTCGCACAGTAACTCCCCAATATGTGAGCATCTGATGTTGCCCGGGCCGAGTTAGTCTTGTGCTCACGGAACTTATTGTATG
>Species-2-samplenameB
AGTAGTGATTTGAAAGAGTTGTCAGTTAGCTCGTTCAGGTAATGGTTCCTCACACTACGTCAAAATAAGAGAGCGGTCGTGACATTATCCGTGATTTTCT
>Species-3-samplenameC
CACTACTATCAGTACTCACGACTCGATTCTGCCGCAGCCACGTATCGCCAGAAAGCCAGTCAGCATTAAGGAGTGCTCTGGGCAGGACAACTCGCATAGT
>Species-3-samplenameD
GAGAGTTACATGTTCGTTGGGCTCTTCCGACACGAACCTCAGTTGGCCTACATCCTACCTGAGGTCTGTGCCCCGGTGGTGAGAAGTGCGCATTTCGTTC

```

I want to split this file in one fasta file per species. I think it’s possible to use the awk function for this but I’m stuck. Does anyone have a script/code that might help me?

Thanks a lot.

# Answer
A simple Biopython solution- iterate over the sequence records, identify the species, open the file handle using append mode to ensure no data is overwritten, and write the record:

```
from Bio import SeqIO

for record in SeqIO.parse("myfile.fa", "fasta"):
    species = record.id.split('_')[0]
    with open(f"{species}.fa", "a") as f:
        SeqIO.write(record, f, "fasta")

```
> 3 votes

---
Tags: python, fasta, bash, awk
---


# Question
Title: Difference between genome assembly and genome sequence alignment to a reference to find structural variants
I'm trying to determine what the difference and benefits of genome assembly and genome sequence alignments are when trying to identify structural variants or transposons in populations. I've been scouring the internet but have only really come across the difference between short vs long reads and de novo assembly vs reference-based.

My understanding is that to identify variations in structural variants within a population there seem to be 2 main comparative genomic methods, the first being what the 1KGP and SDGP did and sequence the whole genome, align the reads to the reference genome and end up with a BAM file.

The second is to assemble personal genomes and then compare or align the assemblies to each other and the reference genome or using the Lastz/LiftOver/ChainNets Examples: 10.1016/j.gene.2005.09.031

Thanks in advance.

# Answer
> 4 votes
> the first being what the 1KGP and SDGP did and sequence the whole genome, align the reads to the reference genome and end up with a BAM file.

If you have a well defined reference genome (e.g. human, mouse etc.) and you are interested in *population level* genetic variation, then this is the main approach. If you sequence a new human genome in the classical way (i.e. short reads ~30X coverage etc.), *de novo* assembly is pointless generally and you can more rapidly map reads to the reference. The read mapping approach has the large advantage that when you get the to variant calling stage, you can use information about both the depth of sequencing and the base quality scores.

> The second is to assemble personal genomes and then compare or align the assemblies to each other and the reference genome or using the Lastz/LiftOver/ChainNets

This is the traditional comparative genomics approach used for non-model organisms and for conducting evolutionary rate comparisons at the *cross-species level*. You discard information on sequencing depth and quality scores, but you can compare many species' genomes at once genome-wide. You also worry less about sequencing errors etc. because you are not looking for relatively rare population level variation, but instead more common (depending on your species) variation between species.

---
Tags: sequence-alignment, assembly, structural-variation, transposable-elements, 1000genomes
---


# Question
Title: Obtaining HGDP project data in fasta format
I need to obtain sample data from modern humans in fasta format. I just need some megabytes of data from every individual. I actually use a script that obtains the cram file from here (ftp.1000genomes.ebi.ac.uk) and then processes it to obtain the fasta file. **The problem is that cram files are large, slow to download and slow to process. It takes days to get the samples. Is there a better way to get these samples in fasta format?**

The script already makes use of samtools to retrieve only the part of the bam file it needs but doesn't help much. Cram files are still gigabytes large for only a few megabytes of data that I need.

I have the same problem with data from the 1000 genomes project.

# Answer
> 0 votes
You can download HGDP data in FASTQ format here: https://www.internationalgenome.org/data-portal/data-collection/hgdp

---
Tags: fasta, samtools, human-genome, 1000genomes
---


title: Convert VCF to genotype table

text:
# Question
Title: Convert VCF to genotype table
How can I convert a VCF file into a genotype table (SNP matrix)?

I have this format:

```
##fileformat=VCFv4.0
##fileDate=20200722
##source=test
##reference=test
##FILTER=<ID=PASS,Description="All filters passed">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
##FORMAT=<ID=DV,Number=1,Type=Integer,Description="Number of high-quality non-reference bases">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  HOR_997 HOR_1020    HOR_1251    HOR_2970    HOR_12367   HOR_473 HOR_8148    HOR_17169   HOR_21258   HOR_4606    HOR_15922   HOR_20226   HOR_9724    HOR_10742   HOR_1071    HOR_10702   HOR_18382   HOR_8819    HOR_9558    HOR_7428    HOR_3023    HOR_18300   HOR_13924   HOR_1080    HOR_4030    HOR_6962    HOR_3460    HOR_8117    HOR_337 HOR_16112   HOR_1626    HOR_8659    HOR_16569   HOR_11083   HOR_18087   HOR_4273    HOR_2237    HOR_8440    HOR_7056    HOR_869 HOR_13987   HOR_12184   HOR_13800   HOR_9808    HOR_10192   HOR_4264    HOR_12830   HOR_10764   HOR_12171   HOR_10784   HOR_610 HOR_11875   HOR_2637    HOR_4136    HOR_282 HOR_10725   HOR_17649   HOR_1384    HOR_2598    HOR_15841   HOR_14342   HOR_18917   HOR_3125    HOR_11662   HOR_6220    HOR_18201   HOR_12412   HOR_14876   HOR_19952   HOR_16358   HOR_6892    HOR_3155    HOR_19184   HOR_9927    HOR_14787   HOR_8370    HOR_4258    HOR_14411   HOR_21548   HOR_12070   HOR_7251    HOR_8904    HOR_4653    HOR_2749    HOR_2957    HOR_2028    HOR_1816    HOR_7236    HOR_495 HOR_2403    HOR_15860   HOR_14439   HOR_11409   HOR_9043    HOR_5876    HOR_13836   HOR_10886   HOR_13938   HOR_20269   HOR_11448   HOR_13979   HOR_14557   HOR_12311   HOR_10004   HOR_15908   HOR_1531    HOR_13821   HOR_18681   HOR_6936    HOR_1990    HOR_10990   HOR_3187    HOR_3728    HOR_15857   HOR_20909   HOR_1581    HOR_13462   HOR_14914   HOR_16149   HOR_12168   HOR_6952    HOR_1048    HOR_18824   HOR_1361    HOR_16071   HOR_3233    HOR_21599   HOR_15871   HOR_11101   HOR_13724   HOR_1780    HOR_8711    HOR_3369    HOR_3909    HOR_21461   HOR_3093    HOR_11013   HOR_18795   HOR_15254   HOR_5486    HOR_17591   HOR_19568   HOR_7129    HOR_6934    HOR_8760    HOR_18924   HOR_21687   HOR_15895   HOR_15898   HOR_15905   HOR_21641   HOR_3686    HOR_7657    HOR_18385   HOR_8710    HOR_17077   HOR_2293    HOR_3081    HOR_7172    HOR_4201    HOR_2383    HOR_19379   HOR_7124    HOR_3912    HOR_15401   HOR_18381   HOR_4034    HOR_16078   HOR_12072   HOR_13716   HOR_4403    HOR_4829    HOR_11431   HOR_4138    HOR_2290    HOR_15788   HOR_7552    HOR_10775   HOR_3153    HOR_14689   HOR_11126   HOR_2180    HOR_2369    HOR_1044    HOR_21428   HOR_18220   HOR_3359    HOR_15778   HOR_15840   HOR_18364   HOR_9250    HOR_11922   HOR_3926    HOR_7474    HOR_12830   HOR_5020    HOR_11450   HOR_9492    HOR_5373    HOR_18913   HOR_12791
1H  42240   .   C   A   999 .   .   GT:DP:DV:GQ 0/0:15:0:100    0/1:9:3:100 0/0:3:0:100 0/0:20:0:100    0/0:4:0:100 1/1:3:3:100 ./.:.:.:.   0/0:24:0:100    0/0:9:0:100 0/0:12:0:100    0/0:14:0:100    ./.:.:.:.   0/0:9:0:100 ./.:.:.:.   0/0:29:0:100    0/0:10:0:100    0/0:5:0:100 0/0:14:0:100    0/0:3:0:100 0/0:18:0:100    0/0:18:0:100    0/0:7:0:100 0/0:9:0:100 0/0:9:0:100 1/1:5:5:100 0/0:13:0:100    ./.:.:.:.   0/0:22:0:100    ./.:.:.:.   ./.:.:.:.   0/0:4:0:100 0/0:18:0:100    ./.:.:.:.   0/0:7:0:100 0/0:11:0:100    0/0:32:0:100    0/0:13:0:100    ./.:.:.:.   0/0:15:0:100    0/0:12:0:100    0/0:8:0:100 ./.:.:.:.   0/0:20:0:100    0/0:11:0:100    0/0:6:0:100 0/0:14:0:100    0/0:20:0:100    0/0:5:0:100 ./.:.:.:.   0/0:3:0:100 0/0:5:0:100 0/0:16:0:100    1/1:3:3:100 0/0:32:0:100    0/0:15:0:100    0/0:17:0:100    0/0:16:0:100    0/0:4:0:100 0/0:17:0:100    0/0:2:0:100 0/0:4:0:100 0/0:7:0:100 0/0:20:0:100    ./.:.:.:.   0/0:20:0:100    0/0:8:0:100 0/0:20:0:100    0/0:3:0:100 0/0:21:0:100    0/0:5:0:100 ./.:.:.:.   0/0:12:0:100    0/0:16:0:100    0/0:87:0:100    0/0:21:0:100    0/0:5:0:100 0/0:20:0:100    1/1:6:6:100 0/0:6:0:100 0/0:10:0:100    ./.:.:.:.   0/0:6:0:100 0/0:4:0:100 0/0:8:0:100 0/0:19:0:100    0/0:17:0:100    0/0:4:0:100 0/0:3:0:100 1/1:7:7:100 0/0:18:0:100    0/0:12:0:100    0/0:11:0:100    0/0:13:0:100    ./.:.:.:.   0/0:18:0:100    0/0:18:0:100    0/0:2:0:100 0/0:12:0:100    0/0:10:0:100    0/0:3:0:100 ./.:.:.:.   0/0:12:0:100    0/0:2:0:100 0/0:8:0:100 0/0:15:0:100    0/0:7:0:100 0/0:22:0:100    0/0:12:0:100    0/0:14:0:100    0/0:8:0:100 0/0:10:0:100    0/0:22:0:100    0/0:10:0:100    0/0:11:0:100    0/0:3:0:100 0/0:6:0:100 0/0:6:0:100 0/0:8:0:100 0/0:18:0:100    ./.:.:.:.   0/0:13:0:100    0/0:20:0:100    0/0:17:0:100    0/0:8:0:100 0/0:4:0:100 0/0:11:0:100    0/0:12:0:100    0/0:8:0:100 0/0:12:0:100    ./.:.:.:.   0/0:6:0:100 ./.:.:.:.   0/0:10:0:100    0/0:6:0:100 0/0:10:0:100    0/0:15:0:100    0/0:15:0:100    0/0:17:0:100    0/0:4:0:100 0/0:9:0:100 ./.:.:.:.   0/0:14:0:100    0/0:13:0:100    0/0:9:0:100 0/0:3:0:100 0/0:12:0:100    0/0:12:0:100    ./.:.:.:.   ./.:.:.:.   0/1:5:2:100 0/0:21:0:100    0/0:9:0:100 0/0:14:0:100    0/0:18:0:100    0/0:11:0:100    1/1:3:3:100 0/0:19:0:100    0/0:25:0:100    ./.:.:.:.   ./.:.:.:.   0/0:19:0:100    0/0:12:0:100    0/0:7:0:100 0/0:10:0:100    0/0:13:0:100    0/0:11:0:100    0/0:13:0:100    0/0:13:0:100    0/0:6:0:100 ./.:.:.:.   0/0:4:0:100 0/0:10:0:100    0/0:9:0:100 0/0:14:2:100    1/1:2:2:100 0/0:16:0:100    ./.:.:.:.   0/0:27:0:100    0/0:18:0:100    0/0:6:0:100 0/0:10:0:100    0/0:14:0:100    0/0:7:0:100 1/1:3:3:100 0/0:15:0:100    0/0:4:0:100 0/0:15:0:100    ./.:.:.:.   ./.:.:.:.   0/0:8:0:100 0/0:16:0:100    0/0:4:0:100 ./.:.:.:.   0/0:12:0:100    0/0:12:0:100    0/0:10:0:100    0/0:7:0:100 0/0:29:0:100    0/0:7:0:100 0/0:8:0:100 0/0:32:0:100
1H  42247   .   C   A   999 .   .   GT:DP:DV:GQ 0/0:15:0:100    0/0:11:0:100    1/1:2:2:100 0/0:19:0:100    0/0:5:0:100 ./.:.:.:.   ./.:.:.:.   0/0:23:0:100    0/0:12:0:100    1/1:13:13:100   0/0:14:0:100    ./.:.:.:.   0/0:9:0:100 ./.:.:.:.   1/1:28:27:100   0/0:12:0:100    1/1:6:6:100 0/0:15:0:100    1/1:3:3:100 0/0:17:0:100    0/0:19:0:100    0/0:7:0:100 0/0:9:0:100 0/0:10:0:100    0/0:4:0:100 1/1:12:12:100   ./.:.:.:.   0/0:22:0:100    ./.:.:.:.   ./.:.:.:.   0/0:6:0:100 0/0:17:0:100    ./.:.:.:.   0/0:8:0:100 0/0:11:0:100    0/0:33:0:100    0/0:14:0:100    ./.:.:.:.   0/0:15:0:100    0/0:11:0:100    0/0:8:0:100 ./.:.:.:.   0/0:21:0:100    0/0:11:0:100    0/0:6:0:100 1/1:15:15:100   0/0:20:0:100    0/0:5:0:100 ./.:.:.:.   0/0:3:0:100 0/0:6:0:100 0/0:17:0:100    0/0:4:0:100 0/0:32:0:100    1/1:15:15:100   0/0:16:0:100    0/0:18:0:100    0/0:5:0:100 0/0:17:0:100    0/0:2:0:100 0/0:5:0:100 0/0:7:0:100 0/0:20:0:100    ./.:.:.:.   0/0:21:0:100    0/0:8:0:100 0/0:20:0:100    0/0:3:0:100 1/1:21:21:100   0/0:6:0:100 0/0:2:0:100 0/0:11:0:100    0/0:15:0:100    0/0:88:0:100    0/0:21:0:100    1/1:5:5:100 0/0:22:0:100    0/0:8:0:100 0/0:6:0:100 0/0:10:0:100    ./.:.:.:.   0/0:5:0:100 1/1:4:4:100 0/0:8:0:100 0/0:19:0:100    0/0:17:0:100    1/1:4:4:100 1/1:3:3:100 0/0:5:0:100 0/0:18:0:100    0/0:11:0:100    0/0:11:0:100    0/0:15:0:100    ./.:.:.:.   0/0:18:0:100    0/0:18:0:100    0/0:2:0:100 0/0:15:0:100    0/0:11:0:100    1/1:3:3:100 ./.:.:.:.   0/0:14:0:100    ./.:.:.:.   1/1:7:7:100 0/0:16:0:100    0/0:7:0:100 0/0:23:0:100    0/0:14:0:100    0/0:14:0:100    0/0:8:0:100 0/0:10:0:100    0/0:24:0:100    0/0:10:0:100    0/0:11:0:100    0/0:3:0:100 0/0:5:0:100 0/0:7:0:100 1/1:8:8:100 0/0:17:0:100    ./.:.:.:.   0/0:14:0:100    0/0:20:0:100    0/0:17:0:100    1/1:8:8:100 0/0:4:0:100 0/0:11:0:100    0/0:13:0:100    0/0:8:0:100 0/0:15:0:100    ./.:.:.:.   0/0:6:0:100 ./.:.:.:.   0/0:11:0:100    0/0:7:0:100 0/0:10:0:100    0/0:17:0:100    0/0:18:0:100    0/0:16:0:100    0/0:4:0:100 0/0:9:0:100 ./.:.:.:.   0/0:14:0:100    0/0:13:1:100    0/0:9:0:100 0/0:3:0:100 0/0:12:0:100    0/0:12:0:100    ./.:.:.:.   0/0:3:0:100 0/0:4:0:100 0/0:21:0:100    0/0:9:0:100 0/0:16:0:100    0/0:16:0:100    0/0:12:0:100    0/0:4:0:100 0/0:20:0:100    0/0:24:0:100    ./.:.:.:.   ./.:.:.:.   0/0:18:0:100    0/0:12:0:100    0/0:7:0:100 0/0:10:0:100    0/0:12:0:100    0/0:12:0:100    0/0:14:0:100    0/0:13:0:100    1/1:6:6:100 ./.:.:.:.   0/0:5:0:100 0/0:13:0:100    0/0:8:0:100 0/0:13:0:100    0/0:4:0:100 0/0:15:0:100    ./.:.:.:.   0/0:29:0:100    0/0:19:0:100    0/0:6:0:100 0/0:9:0:100 0/0:15:0:100    1/1:7:7:100 0/0:4:0:100 1/1:14:14:100   1/1:4:4:100 0/0:16:0:100    ./.:.:.:.   ./.:.:.:.   0/0:8:0:100 0/0:15:0:100    1/1:4:4:100 ./.:.:.:.   1/1:13:13:100   0/0:12:0:100    0/0:10:0:100    1/1:6:6:100 0/0:30:0:100    0/0:8:0:100 0/0:7:0:100 0/0:31:0:100
1H  42264   .   G   T   999 .   .   GT:DP:DV:GQ 0/0:23:0:100    0/0:17:0:100    0/0:3:0:100 0/0:21:0:100    0/0:6:0:100 0/0:8:0:100 ./.:.:.:.   0/0:29:0:100    0/0:20:0:100    0/0:9:0:100 0/0:20:0:100    ./.:.:.:.   0/0:11:0:100    ./.:.:.:.   0/0:26:0:100    0/0:13:0:100    0/0:8:0:100 0/0:19:0:100    0/0:3:0:100 0/0:19:0:100    0/0:18:0:100    0/0:7:0:100 0/0:9:0:100 0/0:18:0:100    0/0:9:0:100 0/0:13:0:100    ./.:.:.:.   0/0:22:0:100    0/0:6:0:100 ./.:.:.:.   0/0:7:0:100 0/0:22:0:100    ./.:.:.:.   0/0:10:0:100    0/0:14:0:100    0/0:39:0:100    0/0:16:0:100    ./.:.:.:.   0/0:21:0:100    0/0:19:0:100    0/0:8:0:100 ./.:.:.:.   0/0:25:0:100    0/0:13:0:100    0/0:6:0:100 0/0:19:0:100    0/0:19:0:100    0/0:4:0:100 ./.:.:.:.   0/0:2:0:100 0/0:11:0:100    0/0:28:0:100    0/0:11:0:100    0/0:37:0:100    0/0:12:0:100    0/0:25:0:100    0/0:19:0:100    0/0:9:0:100 0/0:17:0:100    ./.:.:.:.   0/0:6:0:100 0/0:9:0:100 0/0:28:0:100    0/0:3:0:100 0/0:24:0:100    0/0:11:0:100    0/0:23:0:100    0/0:4:0:100 0/0:16:0:100    0/0:7:0:100 0/0:7:0:100 0/0:19:0:100    0/0:17:0:100    0/0:105:0:100   0/0:23:0:100    0/0:6:0:100 0/0:28:0:100    0/0:18:0:100    0/0:9:0:100 0/0:15:0:100    ./.:.:.:.   0/0:7:0:100 0/0:6:0:100 0/0:11:0:100    0/0:23:0:100    0/0:19:0:100    0/0:7:0:100 0/0:7:0:100 0/0:32:0:100    0/0:20:0:100    0/0:13:0:100    0/0:14:0:100    0/0:18:0:100    0/0:11:0:100    0/0:23:0:100    0/0:19:0:100    0/0:2:0:100 0/0:20:0:100    0/0:15:0:100    0/0:4:0:100 ./.:.:.:.   0/0:17:0:100    0/0:2:0:100 0/0:10:0:100    0/0:17:0:100    0/0:8:0:100 0/0:23:0:100    0/0:20:0:100    0/0:17:0:100    0/0:8:0:100 0/0:14:0:100    0/0:30:0:100    0/0:12:0:100    0/0:15:0:100    0/0:3:0:100 0/0:6:0:100 0/0:9:0:100 0/0:13:0:100    0/0:18:0:100    ./.:.:.:.   0/0:18:0:100    0/0:21:0:100    0/0:23:0:100    0/0:9:0:100 0/0:6:0:100 0/0:14:0:100    0/0:14:0:100    0/0:9:0:100 0/0:18:0:100    ./.:.:.:.   0/0:6:0:100 ./.:.:.:.   0/0:18:0:100    0/0:8:0:100 0/0:15:0:100    0/0:17:0:100    0/0:23:0:100    0/0:26:0:100    0/0:10:0:100    0/0:12:0:100    ./.:.:.:.   0/0:17:0:100    0/0:15:0:100    0/0:11:0:100    0/0:4:0:100 0/0:16:0:100    0/0:15:0:100    ./.:.:.:.   0/0:10:0:100    0/0:7:0:100 0/0:22:0:100    0/0:10:0:100    0/0:21:0:100    0/0:24:0:100    0/0:13:0:100    0/0:19:0:100    0/0:24:0:100    0/0:32:0:100    ./.:.:.:.   0/0:6:0:100 0/0:23:0:100    0/0:17:0:100    0/0:15:0:100    0/0:15:0:100    0/0:15:0:100    0/0:14:0:100    0/0:16:0:100    0/0:17:0:100    0/0:7:0:100 ./.:.:.:.   0/0:8:0:100 0/0:14:0:100    0/0:12:0:100    0/0:30:0:100    0/0:19:0:100    0/0:24:0:100    ./.:.:.:.   0/0:36:0:100    0/0:26:0:100    0/0:7:0:100 0/0:10:0:100    0/0:21:0:100    0/0:8:0:100 0/0:13:0:100    0/0:17:0:100    0/0:4:0:100 0/0:18:0:100    0/0:2:0:100 ./.:.:.:.   0/0:9:0:100 0/0:18:0:100    0/0:5:0:100 ./.:.:.:.   0/0:9:0:100 0/0:14:0:100    0/0:15:0:100    0/0:11:0:100    0/0:32:0:100    0/0:9:0:100 0/0:12:0:100    0/0:34:0:100

```

And I need a SNP matrix in this format:

```
 , pos1, pos2, pos3
ID1   0,    1,    0
ID2   1,    1,    0

```

# Answer
If you are familiar with plink and python, you can use pandas\_plink.

```
import pandas_plink
snp_info,sample_info,genotypes  = pandas_plink.read_plink('genotypes/chr.1')
genotype_mat = genotypes.compute()

```

`genotype_mat` is a matrix of genotypes (in 0,1,2) and `sample_info` has the sample names and `snp_info` has the SNP positions, ref, alt etc. You need to convert vcf to plink format first. Please make sure ref and alt alleles are extracted as expected.
> 1 votes

---
Tags: python, vcf, snp, gwas
---


# Question
Title: Download or create VCF from Human Genome Diversity Project (HGDP)
I can't seem to find HGDP data in VCF format.

Does anyone know where I can download it from?

Perhaps there a Snakemake pipeline for calling VCFs from this data?

# Answer
I presume you are looking for the recently sequenced WGS data rather than the genotype data? If so, the different files can be found here: ftp://ngs.sanger.ac.uk/production/hgdp/hgdp\_wgs.20190516/

Each chromsome has a seperate vcf file with a corresponding .tbi index file.

The gvcfs can be found here: ftp://ngs.sanger.ac.uk/production/hgdp/hgdp\_wgs.20190516/gVCFs/

These files are more compressed and smaller in size than the regular vcfs.
> 3 votes

---
Tags: vcf, variant-calling, data-download, human-genome
---


# Question
Title: europmc and tidypmc R library for extracting or making metadata from publication
I'm trying to use this europmc r libray where I have a list of pmids to look for. I tried with pubtator but its bit complicated.In Europmc i can all the annotated terms etc.

```
library("europepmc")

example list of PMIDS
30024784
30555165
30510081
31688884
31516032
28588019
29286103

```

Now what I'm doing is looking each ID using epmc\_details function which is not the way i would do if i have to look for hundreds.

```
epmc_details(ext_id = '30510081')

```

I have to question how can I run the epmc\_details in a loop or some other way where i can look for PMIDS one by one and save the result in a data frame.

The epmc\_details is returned as a list. The structure of the list is as such

```
[1] "basic"           "author_details"  "journal_info"    "ftx"             "chemical"        "mesh_topic"      "mesh_qualifiers"
[8] "comments"        "grants"

```

I would only like to save basic,chemical,mesh\_topic,mesh\_qualifiers in one data frame.

For example if my first id is this `30510081` the dataframe should have first column as my ID which is basically `basic[1]`and rest of the information appended to the next columns.

such as

```
 ID chemical mesh_qualifiers mesh_topic gene

```

Any suggestion or help would be really and highly appreciated

I was looking at europmc site through the browser

this was one of my query

when i highlight the keyterms i do see in the abstract itself all the keyterms are getting annotated but when i do the same query search through R I do see empty results as such why there is a difference?

```
$chemical
# A tibble: 0 x 0

$mesh_topic
# A tibble: 0 x 0

$mesh_qualifiers
# A tibble: 0 x 0

```

I found better way of getting data from pubmed using the tidypmc library.

```
library(tidypmc)
doc <- pmc_xml("PMC6365492")
doc
txt <- pmc_text(doc)
txt

count(txt, "section")

cap1 <- pmc_caption(doc)
filter(cap1, sentence == 1)

tab1 <- pmc_table(doc)

sapply(tab1, nrow)

tab1[[1]]
attributes(tab1[[2]])

collapse_rows(tab1, na.string="-")

library(tibble)
x <- xml_name(xml_find_all(doc, "//*"))
tibble(tag=x) %>% count("tag")

library(tidytext)
x1 <- unnest_tokens(txt, word, text) %>%
  anti_join(stop_words) %>%
  filter(!word %in% 1:100)
#  Joining, by = "word"
#filter(x1, str_detect(section, "Case description"))

filter(x1, str_detect(section, "Results"))
count(a$word)

tbls <- pmc_table(doc)
map_int(tbls, nrow)

tbls[[1]]

collapse_rows(tbls, na.string="-")

```

But if i understand it can use one PMC id at a time. Again keeping my original question how can i put this in a loop to query lets say i have 100 PMCID and get it result and store in a dataframe.

After using tidypmc i found that i parse all the publication based on attributes or tags. Such as title,abstract, results etc etc.

Lets say I'm interested in the table tags information where they have metadata of patients as well as others.So if a paper contain multiple tables I would like to store each of them in a data-frame under the respective publication.Since I have multiple IDs to search and do the save as above mentioned. How to put this through a loop or can it be done without using loop ?

Any suggestion or help would be really appreciated as always.

# Answer
> 3 votes
The idea for this code is to first convert PIDs to PMCIDs, then run the `tidypmc` in a loop over the PMCIDs. The only problem is that `tidypmc` failed to retrieve tables from most of the IDs in your example list.

```
library(tidyverse)
library(tidypmc)
library(httr)
library(jsonlite)
example_pids <- c(30024784, 30555165, 30510081, 31688884, 31516032, 28588019, 29286103) %>% as.character()

#-- Convert to PMC ids
convertPIDtoPMCID <- function(pids) {
    #-------- Make API request
    pids4query <- paste(pids, collapse = "%0D%0A")
    idconv_req <- paste0("https://www.ncbi.nlm.nih.gov/pmc/utils/idconv/v1.0/?ids=", pids4query, "&idtype=pmid&format=json&versions=no&showaiid=no&tool=&email=&.submit=Submit")
    pids_json <- GET(idconv_req) %>% content("text") %>% fromJSON()
    #-------- Get info from JSON
    pid2pmc <- pids_json$records %>%
        select(pmcid, pmid) %>%
        as.data.frame()
    rownames(pid2pmc) <- pid2pmc$pmid
    pmcids <- pmcids <- pid2pmc[pids, "pmcid"]
    return(pmcids)
}
example_pmcids <- convertPIDtoPMCID(example_pids)

#-- Try to get data with tidypmc
pub_tables <- lapply(example_pmcids, function(pmc_id) {
    message("-- Trying ", pmc_id, "...")
    doc <- tryCatch(pmc_xml(pmc_id),
                    error = function(e) {
                        message("------ Failed to recover PMCID")
                        return(NULL)
                        })
    if(!is.null(doc)) {
        #-- If succeed, try to get table
        tables <- pmc_table(doc)
        if(!is.null(tables)) {
            #-- If succeed, try to get table name
            table_caps <- pmc_caption(doc) %>%
                filter(tag == "table")
            names(tables) <- paste(table_caps$label, table_caps$text, sep = " - ")
        }
        return(tables)
    } else {
        #-- If fail, return NA
        return(NA)
    }
})
names(pub_tables) <- example_pids

#-- Inspect results
pub_tables$`30555165`$`Table 1 - Patient demographic and baseline characteristics`
pub_tables$`29286103`$`Table I - Sample summary.`

```

Tables will require quite a bit of tidying after this. Good luck!

---
Tags: r, public-databases, parsing
---


# Question
Title: How to retrieve the DNA sequence from a particular species and region through command line?
I am new to working with the NCBI database so I am not sure how trivial this question is. I wanted to get DNA sequence of a specific version of the Xenopus tropicalis genome of a specific region (ie 1-100bp) without doing it manually each time (going on genome browser). I know there is a faster way and efficient way to do this through the command line but I found the NCBI tutorial hard to follow.

# Answer
**Samtools**

Samtools is a toolkit that includes a binary called `samtoools`, which can be used to query a FASTA file — to extract a subsequence based on some coordinate range you specify.

If you do not have `samtools` installed, you need to add it to your computer. For OS X, you could use Homebrew, for example:

```
$ brew install samtools

```

For Ubuntu, you could use `apt-get`:

```
$ sudo apt-get install samtools

```

**Reference assembly**

The next step is to download a reference assembly for your genome and genomic build-of-interest.

There is a research group dedicated to making *Xenopus* datasets available. If you were interested in v10 of the `Xentr` genome, for instance, you could do the following or similar to download it:

```
$ curl --output "XENTR_10.0_genome.fasta.gz" "http://ftp.xenbase.org/pub/Genomics/JGI/Xentr10.0/XENTR_10.0_genome.fasta.gz"

```

You could also dig through the NCBI site to get a similar link. It's the same idea — just change the web address, accordingly.

**Indexing**

Indexing makes it possible to extract an arbitrary subsequence of interest from this FASTA file. You need to index, but you only need to do it once:

```
$ samtools faidx XENTR_10.0_genome.fasta.gz

```

This creates two new files: `XENTR_10.0_genome.fasta.gz.fai` and `XENTR_10.0_genome.fasta.gz.gzi`.

*Note*: You need to keep these two files in the same directory as the original, compressed FASTA file `XENTR_10.0_genome.fasta.gz`.

**Querying**

Now you are ready to query, or extract a subsequence of interest. If you wanted the first hundred bases of the first chromosome, for instance:

```
$ samtools faidx XENTR_10.0_genome.fasta.gz Chr1:1-100
>Chr1:1-100
cctcagcgtcttcttctttctgctcccactgcattagcatagggagagagggccacagag
gctaggtagcattgccggctagcaatcttcagcgtcctct

```

The search term used here is `Chr1:1-100`. Replace `Chr1` with the chromosome of interest and `1` and `100` with the start and stop coordinates, respectively, of interest. Coordinates used with `samtools` queries use a one-based indexing scheme.

*Note*: This uses the chromosome naming scheme in the `XENTR_10.0_genome.fasta.gz` file (`Chr1`, etc.). Other research groups (NCBI) might use different chromosome names. If you get your FASTA from somewhere else, investigate the FASTA file to determine the naming scheme being used.

If you want to redirect output to a file, use the output redirect operator (`>`):

```
$ samtools faidx XENTR_10.0_genome.fasta.gz Chr1:1-100 > mySubsequence.fa

```
> 4 votes

# Answer
I'd download the fasta, and use samtools, or maybe bedtools.
> -2 votes

---
Tags: phylogenetics
---


# Question
Title: Assign cell types to groups of cells based on their gene expression profiles
I have large filtered, normalized dataset of `scRNA-seq` data of `C.Elegans` species. Rows are genes (10 000), columns are cells (66 000). Let's say that I got 40 different groups of the cells based on their expression profiles, how could I now assign cell types to these groups? I would expect to get some reference dataset for `C.Elegans` and compare it with what I got, but I honestly have no idea of where to get such a dataset, and whether there are any `R` packages that do that. Any advice would be greatly appreciated.

# Answer
> 4 votes
I think the analysis will involve two steps:

**1) Compile a reference list** of cell types, with a list of genes expressed in each of them, based on the literature. *C. elegans* has only ~1000 cells, so I'm sure it's been well researched. Mark which genes are just expressed in that particular cell type, and which are uniquely expressed, this latter will be more useful markers.

Resources to do that:

This will involve mostly manual work, but maybe you can find some list.

For example, platelets express *VWF*. And activated platelets express *SELP*.

**2) Use one of the single-cell RNA-seq packages** (e.g. Seurat) to cluster cells and to find the "marker" (DE) genes for each cluster (vs. all other clusters).

For each cluster, look up the marker genes in the reference list and infer the cell type.

You may need to infer the cell type from the functions and descriptions of the genes. Also see Retrieve detailed gene descriptions

Keep in mind that you will have many confusing doublets in the dataset expressing characteristics (genes) of two cells.

You may also need to re-do t-SNE and clustering on a subset of cells, e.g. blood/immune cells, in an iterative, hierarchical way.

---

Do you have information on which cells are from the same worm? Once you have a reference list, you can try and identify cells of just one worm.

# Answer
> 0 votes
This is an old post, but we have recently published a web-based tool called CIPR (cluster identity predictor) that enables users to upload their data along with a custom reference file for annotating unknown clusters. You can probably obtain relevant reference datasets from GEO.

You can read the manuscript in BMC Bioinformatics for more information about how the algorithm works.

---
Tags: r, scrnaseq
---


# Question
Title: Specific cell type identification in Single Cell Sequencing
In order to define which cell is of which type we need to identify a set of rules, for instance neurons should express one of the following: `Thy1`, `Rbfox3`, `MAP2`, `Camk2b`, `Gad1`,`Cck`, `Reln`, and should not express any of the following: `cd45`, `Tmem119`, `CD11b`, ... and others. Maybe `Rbfox3` should always be expressed in any neurons. So, these rules need to be compiled manually as I understand to identify each cell type, from the literature. I would assume that many people already faced this issue, and maybe came up with software to do that? Are there any software where we would be able to supply a matrix of single cell expression data and have the types assigned with certain degree of certainty?

I know that maybe we can also use heatmaps for that, like:

cell type specific genes - heatmap using rank based approach

But again, this seems too manual work. Besides, I do not think that we can assign different neuronal types based on just a heatmap.

# Answer
> 6 votes
I do not know of such software.

However, I believe this effort is a bit misdirected. The purpose of single-cell sequencing is to get a better understanding of cells; their heterogeneity and functional diversity or developmental / biological processes such as differentiation, using a higher "resolution" method. In other words, if we had the methods you are asking for, there would be no need for the single-cell experiment.

Running your data through some pipeline utilizing previous knowledge has the danger of forcing ideas on the data, rather than seeing what the data tells you: it would be better to try and understand and explain the biologically relevant heterogeneity and diversity in your cells, together with a critical comparison of cell type characteristics (expression of genes, pathways) with knowledge described in literature. This question is also very similar to the one you previously asked and was answered.

Also good to recognize that your approach has a lot of implicit assumptions, such as: cells can be clearly categorized as +/- expression for each gene, or it ignores the systems level (networks and pathways), thus again forcing concepts on the data.

---

This said, one method would be to construct a table listing cell type markers as described in the linked answer above, then write a script that determines a cutoff value for expression of genes in your data (see this), then ranks cell types for each cell. For example, you can measure the number or proportion of uniquely genes expressed (although that's very simplistic).

Another option, if you are familiar with machine learning, is to train a classifier on an annotated dataset, and then use that on new data.

Also see a convenience function below, which requires a table of marker genes for each cell cluster (i.e. (the `Seurat::FindAllMarkers()` output), and a reference df listing genes and corresponding cell types in its HGNC\_symbol and Cell\_type columns, and returns the table with listing clusters, their marker genes and corresponding cell types.

```
getCelltypes <- function(markers, reference) {
  marker.celltype <- markers
  marker.celltype$Cell_type <- marker.celltype$gene
  marker.celltype$Cell_type <- with(reference, celltype[match(marker.celltype$Cell_type, gene)])
  return(marker.celltype)
}

```

Many variants of this reference table and function can be created, feel free to modify it.

In the future, I think we can expect such classifier software that uses reference data from the Human Protein Atlas, HCA and similar projects.

# Answer
> 1 votes
SingleR, a novel computational method for unbiased cell type recognition of scRNA-seq: https://github.com/dviraran/SingleR

# Answer
> 1 votes
**singleCellNet** is a computational method using supervised machine learning for **quantitatively** **cell type annotation**. singleCellNet also enables **cross-platform** and **cross-species** comparison. It is also actively-supported.

Here is the github https://github.com/pcahan1/singleCellNet , and the vignette https://pcahan1.github.io/singleCellNet/

# Answer
> 0 votes
We have recently published a web-based tool called CIPR (cluster identity predictor) that helps with cluster annotations. You can select one of 7 reference datasets (some contains neurons) or upload a custom reference gene expression file for predicting cluster identity. You can read the manuscript for more information about how the algorithm works, but in a nutshell, it compares differentially expressed genes or global expression patterns within your data to known reference datasets and creates quick visual outputs. I find it helpful for interactive analyses where I try different clustering parameters and such.

This article is a good resource that compares different cluster calling algorithms.

---
Tags: scrnaseq
---


# Question
Title: Is there a computational method/method/way to predict if a cell population in a tissue immigrated or is enriched locally?
I have single cell data which I have analysed for differential expression. In my experiment, I subjected two groups of mice (a control and a treatment group) to treatment that will lead to the proliferation of my cells of interest in the treated mice. Other scientists have said that the increase of cell numbers caused by this treatment results from increased recruitment from precursor cells from the bone marrow while other scientists argue that the increase in cell number is due to local proliferation.

I wish to address this question computationally, based on my single cell rna seq data using R or any other method. I have googled but I can find no suitable tool. Does it suffice to test the expression of known cell migration markers, produce a migration score and judge based on the migration scores of my proliferating cells?

If not, I will be very grateful for suggestion of a computational tool, method or package that can help me address this question even indirectly.

Thank you in advance.

# Answer
> 1 votes
This may not be a viable answer for all experimental contexts, but you can compare expression profiles of individual cells against datasets that contain known reference cell types (migratory and not). We have recently published a web-based tool called CIPR (cluster identity predictor) that helps with cluster annotations in scRNAseq experiments. There are various reference datasets available, but you can also provide a custom reference by preparing data from databases like GEO.

You can read the manuscript for more information about how the algorithm works, but in a nutshell, it compares differentially expressed genes or global expression patterns within your data to known reference datasets and creates quick visual outputs. If your reference contains both migratory and non-migratory cell types, maybe CIPR can help differentiate them. For other identity prediction tools, check out this article is a good resource that compares different cluster calling algorithms.

---
Tags: scrnaseq, bioconductor, software-recommendation, computation
---


title: no result from heat map WGCNA

text:
# Question
Title: no result from heat map WGCNA
I have created a heat map using WGCNA using R. But the heat map command is not showing the plot and there is no error. I have a large number of modules and traits. I am not sure if that is the reason. But I tried making it in HPC R and it still is not showing any result plot. Thank you!

# Answer
> 1 votes
The plot is probably being drawn to a device that you cannot view on a cluster. If you're using the `plotNetworkHeatmap` function, say, try something like this:

```
png("plot.png", ...) #specify height, width, etc
plotNetworkHeatmap(....) #call your function
dev.off()

```

Now, `scp` the `plot.png` file over to your local machine to view the plot.

# Answer
> 1 votes
If you've copied over from some of the WGNCA tutorials you might have additionally copied the sizeGrWindow command - it's caused some of my plots to disappear.

ie:

```
sizeGrWindow(9,6)  # delete this line
plot(...)

```

---
Tags: r, hpc
---


# Question
Title: How can I take cell number into account to find total RNA expression?
I'm hoping to quantitatively show differences in *total* RNA expression for a gene in a cluster of interest between different experimental groups. My exported *average* RNA values for each experimental group in my macrophage cluster are shown below:

00 d Control: 10.43

00 d PLX5622: 14.48

07 d Control: 16.08

07 d PLX5622: 17.12

28 d Control: 15.90

28 d PLX6722: 19.40

However, the violin plot for Cd68 in the macrophage cluster looks like this

> VlnPlot(combine.combined, features = ("Cd68"), pt.size = 1, idents = "Macrophages", group.by = "orig.ident", split.by = NULL, assay = "RNA"):

The violin plot suggests that the Control groups actually have more total CD68 RNA in them than in the PLX5622 groups, which isn't what the average RNA values show. Any code you can share showing me how to get total RNA per group would be an amazing help!

Thanks so much!

# Answer
(I currently cant comment as I used up all my rep for a bounty.)

Could you explain how you calculate the average? In you plot the averages are between an expression level of 2 and 3, while you report values over 10. Given the data you show, there is no significant difference and I would guess p-values are above 0.8!

And you mention that the violin plots suggest the control groups have higher expression.. but how do you see that? For me they are identical with minimal trends that the PLX5622 groups are higher.

And to address your question regarding showing the total: No, you would want to show the mean or median, as the total is pointless, given deviating sample sizes.

EDIT: Overall, I have a feeling you uploaded the wrong picture!

EDIT2: I think I got your point now: The difference in mean is caused by an increased number of outliers in your controls at the value 0. Despite those, the means are identical and this is what can be seen through the violin plots. You should rather try to diminish the effect of these outliers.
> 1 votes

# Answer
To me, it looks like the widest parts of the plots are higher in the treated. I don't see why you are so sure that the violin plots are hugely different from the calculated averages.

Either way, the differences are really small. Maybe not statistically significant.

BTW, in the future, you might want to disguise your treatment and gene of interest when on a public forum.
> -1 votes

---
Tags: seurat, single-cell
---


# Question
Title: using snakemake shadow rules to store temp files on local nodes
Snakemake is not really designed to store temporary files to local disk, but in the documentation, they suggest a workaround using "shadow rules".

The problem I have is, that our front node does not have `/scratch/..` directory and therefore snakemake complains right away

```
snakemake my_result.tsv -p --jobs 10 --shadow-prefix /scratch/$USER --cluster "qsub "
   ...
   PermissionError: [Errno 13] Permission denied: '/scratch'

```

Every single computational node does have the directory, only the frontend does not. Is there any easy workaround?

--- edit ---

I will try to clarify where I think is the problem.

The head node does not have anything mounted to `/scratch`, I think that's why I am getting that error. I sort of expected that snakemake will attempt do create the shadow environments for individual executed jobs, because if it does, everything would be fine. All the computing nodes do have their local disks mounted to `/scratch`.

Of course, I could use a different shadow directory, but all the directories that are mounted both at the front-end and computing nodes are NOT local disks of individual nodes, which would cause more network traveling of the temporary data.

--- edit 2 ---

Thanks @Chris\_Rands for finding the relevant lines. I opened an Issue on Snakemake GitHub. Once I will know more, I will get back with the answer here.

# Answer
What is happening is your Admin is not allowing you access to /scratch ... which is undertandable. Essentially you are attempting to access the headnode `/scratch`. If every user on the cluster had access, it would be a bit chaotic.

I am surprised by the error, not that I would use Snakemake I'd code it, because you would expect a data pipeline tool as established as Snakemake to have a solution to working on a cluster. What you are describing is a standard cluster set-up between headnode and 'sub-node' (the technical term is slavenode, but I hate that term). The headnode is mirrored across the cluster, hence the error. Once you submit via `qsub` the mirror will 'display' the headnode layout on the subnode.

* You could ask your admin for permission to `/scratch`, this is root privilege. It is a bit unlikely they will say 'yes' but its worth a try.
* Otherwise what about a virtual environment? When Perl suffers this problem I simply use 'perlbrew', anaconda has very nice virtual environments that are easy to setup
* @Chris\_Rand points out piping to a different output.
* I bet there is a solution in '.profile' and thats how I solved comparable problems in the past.
* Of course `Snakemake` do a Q&A interaction on Overflow SE (but please come back here!) ... they must know the work around.

Frankly I am astonished because the whole point of Snakemake, is the ease of constructing a data pipeline and data pipelines are processed on clusters, thus this error seems to up-end it, but who am I?
> 0 votes

---
Tags: snakemake
---


# Question
Title: How to generate Beta diversity boxplot from phyloseq object?
I am working with profiled metagenomic taxonomy abundance data. I want to generate beta diversity (bray-curtis) boxplot from a phyloseq object where two groups (control and test) will be shown. Something like this figure:

How can I make it? Can anyone share any code or tutorial?

Thanks, dpc

# Answer
> 1 votes
You could convert the phyloseq object to a dataframe and just plot it using base R or ggplot2. I'm not exactly sure what the bray curtis df will be called but you can access it using the @ notation then convert with data.frame e.g.

```
bray_for_plot <- data.frame(phylo_obj@bray)

```

---
Tags: statistics
---


# Question
Title: Converting aligned fasta to plink ped/bed
I have an alignment of multiple sequences in a FASTA file (output from MAFFT), for which I would like to simulate a phenotype using plink, but for that I need to have my alignment in a PED file or plink BED. There's plenty online about converting in the other direction, but the only mention I found of converting FASTA to BED is about aligning a raw FASTQ file for use in plink - my data is already aligned! How can I do this conversion?

Note that How to convert FASTA to BED is about the UCSC BED format, not the plink binary PED.

# Answer
An alignment can be the result of two slightly different analyses:

There is **multiple sequence alignment** (which is what you get from MAFFT) where sequences are aligned so that similar regions are on top of each other. This may require introducing indels (insertions / deletions) if a particular region is absent in some of the sequences.
Then there is **alignment to a reference sequence**. This is usually how you get variant data and most of the tools deal with this use case. Therefore, in the link you provided the solution includes aligning the data to a reference sequence. Usually the refernce sequence is much larger (e.g. the human genome) and after alignment the result is a file (.bam) that tells you where the query sequences match in the reference.

I did a quick search for ways of converting a multiple sequence alignment (MSA) to VCF. There is a tool called `msa2vcf` in the Jvarkit collection of utilities that can do this. The example is for the CLUSTAW format but FASTA is accepted as well.
You shoud check that indels were treated correctly because these are most likely to cause trouble. Then you can simply convert the VCF to PLINK format using PLINK.
> 2 votes

---
Tags: fasta, format-conversion, plink
---


title: Seurat DE t.test

text:
# Question
Title: Seurat DE t.test
I am new to the Seurat package and was looking at the tutorials on the Seurat website and am just curious to know if there is any way to perform a t.test for the expression of a specific gene between two groups within the Seurat package, even possibly showing these results with box plots? Or any way to extract this information from the Seurat object in order to do this using ggplot or anything of that sort? I don't have any specific examples as I've just been going through the Seurat tutorials trying to understand everything, but if anyone has any sort of advice on the topic I would greatly appreciate it!

# Answer
Most of your questions are part of the Vignettes on the Seurat home page.
In the differential expression Vignette you can see that you can run Student's t-tests in the `FindMarkers` wrapper function.
However, simple t-tests may not be the best option for single cell data. There are many algorithms for DE testing with some being directly accessible through `FindMarkers`. The alogrithm implemented in `MAST` is recommended frequently.

Regarding the box plots you should have a look at the violin plot function of Seurat:

```
VlnPlot(SeuratData, c("gene_of_interest"), group.by = "groups_compared_in_DE_test")

```
> 1 votes

---
Tags: gene, proteins, genome, seurat
---


# Question
Title: Seurat FindMarkers() output interpretation
I am using FindMarkers() between 2 groups of cells, my results are listed but i'm having hard time in choosing the right markers. Do I choose according to both the p-values or just one of them? If one of them is good enough, which one should I prefer? I am interested in the marker-genes that are differentiating the groups, so what are the parameters i should look for?

```
             p_val avg_logFC pct.1 pct.2  p_val_adj
UBD   3.080608e-06 1.7755753     1 0.000 0.07683959
CFB   1.262305e-05 1.3989233     1 0.067 0.31485675
RGS13 1.548593e-05 0.9009480     1 0.200 0.38626551
PYDC1 1.896309e-05 0.9622537     1 0.133 0.47299636
BIRC3 4.099414e-05 1.0129472     1 1.000 1.00000000
ICAM1 4.379895e-05 0.8219610     1 0.200 1.00000000

```

as you can see, p-value seems significant, however the adjusted p-value is not. and when i performed the test i got this warning `In wilcox.test.default(x = c(BC03LN_05 = 0.249819542916203, ... : cannot compute exact p-value with ties` I am completely new to this field, and more importantly to mathematics. Please help me understand in an easy way.

# Answer
> 1 votes
The p-values are not very very significant, so the adj. p-value. You need to plot the gene counts and see why it is the case. It could be because they are captured/expressed only in very very few cells. `VlnPlot` or `FeaturePlot` functions should help. You haven't shown the TSNE/UMAP plots of the two clusters, so its hard to comment more.

---
Tags: r, scrnaseq, seurat, p-values
---


# Question
Title: Help with setting DimPlot UMAP output into a 2x3 grid in Seurat
On my merged seurat object of 6 samples, when I use the split.by function in tandem with the Dimplot/UMAP plot, all six samples are displayed in series along a commonly labeled 'UMAP\_1' x-axis in an order of (these are arbitrary for simplicity, here) Z, Y, X, C, B, A. Instead, I would like to specify these plots to output as a 2 row x 3 column arrangement ordered (1st Row) A, B, C and (2nd Row) X, Y, Z.

I'm assuming I could use something like: group\_levels \<- c("A", "B", "C", "X", "Y", "Z") to define the order I'd like, but I'm lost on how to get this to arrange the UMAP plots into a 2x3 grid I'd like.

Appreciate any help!

# Answer
> 2 votes
> how to get this to arrange the UMAP plots into a 2x3 grid

You can use the `ncol` parameter, which will set the number of columns in the "combined plot".

And for future reference, if you get the help page of the function you are using with `?DimPlot`, you will see the `...`. This means that you can feed additional arguments to other functions that `DimPlot()` uses. Apparently,  `CombinePlots()` is used under the hood and has an `ncol` parameter.

---
Tags: seurat, single-cell
---


# Question
Title: Identification of starting point for the synthesis of a protein
# identify where is the starting point for the synthesis of a protein in the DNA

For indentify a starting point in RNA or DNA my function search all ATG or M in the sequece and report the linear position of start, than the function search TAA and TAG for stop codon and report the end position.

For the end codon what I have to search? only TAA and TAG or TGA and TGG too?

How is possible some protein begin with "F" as in "FFRE..." and not "M..." ?

Than it result more starting point solutions respect to protein found in GeoBank, what are the add restriction for limit that number?

Are there easy to search sequences that make sure that is a start of the protein code? as seems to me "AGG AGG" but is not clear where they should be in the array and not seems general seen proteins.

Can the synthesis of protein be affected by temperature? What other factors come into play?

Passing my function to DNA of one virus, the first test, the function identifies in the end of virus many little proteins of 2 3 4 5 10 protein letters. These little protein pieces are really sintetized? or there is a min lenght for sintetized proteins?

# Answer
> 1 votes
You seem to be interested in identifying potenital coding regions in a DNA sequence. These regions are called open reading frames (ORF) because they are open to translation.
As detailed in the linked Wikipedia article, an ORF usually starts with an ATG and ends with one of the stop codons (TAA, TGA or TAG in DNA space). However, there are many ORFs that are not translated at all, particularly short ones.

Since finding ORFs is an old and established task in bioinformatics there are many algorithms availabe to do this. I would suggest using such an existing tool because they often have additional features and bugfixes. There is a list in the link above.

Estimating the stability and likelihood of a particular amino acid sequence is much more complicated. There are many tools for protein structure prediction that may be worth having a look at.

---
Tags: proteins, covid-19, sars-cov-2
---


# Question
Title: How to get transcriptome FASTA file for viruses for Kallisto pseudo-alignment?
I guess the title is self-explanatory. I'd like to find the reference transcriptome files for a few human viruses. My RNA-seq samples are from human tumor cells that were infected with viruses. The main part of my analysis is quantifying human gene expression from these samples. However, I also would like to find what proportion of my raw data is originated from virus RNA.

So, I'd like to learn how to find the reference transcriptome data for human viruses. Moreover, I appreciate it if you could give me some insight into your workflow for approaching this problem. Finally, if there is any specific command line, R, or python tool for this specific task, I'd like to know about it.

# Answer
> 1 votes
You might be interested in the viral sequence data from CATCH. There are various different viral subsets available.

You could also try the NCBI virus website, filtering on human hosts. There's a "Download" link on this page that allows you to download the full sequence set:

https://www.ncbi.nlm.nih.gov/labs/virus/vssi/#/virus?VirusLineage\_ss=Viruses,%20taxid:10239&HostLineage\_ss=humans,%20taxid:9605&SeqType\_s=Nucleotide

You've asked a very specific question, so I'm giving you a very specific answer. It seems that what you *actually* want to do is different from the self-explanatory question.

For the purpose of creating an index for Kallisto or Salmon, it should be okay to use the full dataset. Those programs are designed to handle eukaryotic transcriptomes as input, which typically include many gene copies and isoforms, so transcript / genome duplication shouldn't be much of an issue.

Without knowing more about your particular application, it's difficult to provide additional advice. From what you've written, i.e. "what proportion of my raw data is originated from virus RNA", a hit to *any* viral sequence would be sufficient. For example, if you want any more specific information about viral proportion / species, then you should be using a metagenomic tool like Kraken2, rather than Kallisto / Salmon.

# Answer
> 1 votes
I found a possible answer to my own question by reading this biostar page and also finding an NCBI page for downloading the reference genome file for viruses.

* You should go to this NCBI page and search for your target virus.
* Then click on the virus name and you will be redirected to a page like this which gives you an Entrez records table. You should click on the direct link in the genome row in the table.
* Click on `Genome Assembly and Annotation report`  on the new page. You will be redirected to a page like this with some information on reference files submitted to NCBI. The last column in this table is for links to the NCBI FTP site. There are two options to access RefSeq or GenBank FTPs. You can choose either of them. I chose RefSeq. Then, you will be redirected to a page like this with different fasta files for that specific virus.
* the file with the name `GCF_002815995.1_ASM281599v1_cds_from_genomic.fna.gz` is the transcriptome file annotated with CDS info which can be used for Kallisto alignment.

**Note**: If you couldn't find your desired virus genome in the NCBI page mentioned above, you can proceed with the method for generating your desired transcriptome file from genome and annotation file using `gffread` package. The details on this method are provided in the answer to the Biostar question mentioned above.

---
Tags: rna-seq, covid-19, sars-cov-2, kallisto, virus
---


# Question
Title: Understanding the fasta file format and using edit distance for alignment
I'm a computer scientist teaching algorithms development in the Fall. One of the algorithms we teach is called Edit Distance, and our folklore is that it is used to compare RNA sequences (is this actually true in practice?).

I would like to have students implement the edit distance algorithm and run it on actual SARS-COV-2 sequences, so I'm trying to understand exactly what I get from the GenBank database. I downloaded this one: https://www.ncbi.nlm.nih.gov/nuccore/1798174254

I am looking at the genomic.fna file. So this is apparently the FASTA file format, and the lines that begin with \>MN988669.1 ... are comments. I see comments like:

```
>MN988669.1 Severe acute respiratory syndrome coronavirus 2 isolate 2019-nCoV WHU02, complete genome

```

Followed by an RNA string. Is this the start of a new sequence for a different coronavirus specimen? So I could have students extract each of these and run edit distance and then produce a dendrogram or something? How do I find more information on where the samples came from? Is this the right file to be using, or should I use the gbff file? And are the PDB files interesting to me at all (I actually do know what PDB files are)?

Also, are there any recommended data sets where we could do something like track mutations in the virus (and see, e.g., that the NYC outbreak originated from Europe and not China)? Are there other useful algorithms / data that might be interesting for students to study in this vein? Specifically interesting to me would be graph-search algorithms, minimum spanning trees, and network flow. Also any NP-complete algorithms that we could run backtracking on. Obviously taking the theoretical study of algorithms to something as currently topical as the coronavirus has pedagogical value.

Thanks

EDIT:

Based on comments below here is what is taking shape.

1. Have students implement vanilla EditDistance (which there seems to be some disagreement about which algorithms are named what, so let's say insertions and deletions only, which I would call Longest Common Subsequence LCS). Then a variant that also does alignment (i.e. full Levenshtein distance computation, which I would call EditDistance, but Wikipedia calls the Needleman-Wunsch algorithm with a gap penalty of 1), then maybe Needleman-Wunsch with different gap penalties (if someone tells me what would make sense biologically).
2. Have students implement basic heirarchical clustering / phylogenetic tree generation algorithms *a la* https://www.ncss.com/wp-content/themes/ncss/pdf/Procedures/NCSS/Hierarchical\_Clustering-Dendrograms.pdf.
3. Have students run their sequence alignment variants and different clustering algorithms on SARS-COV-2 sequences and report on how the choices of parameters in 1 and 2 change the result and therefore potentially the analysis.
4. Ask some open ended written response questions on what this might mean for society, whether this introduces ethical considerations for an algorithms designer or are they doing just the math, etc.

My learning objectives (as they are now shaping up) are:

1. Students will understand that just because their algorithm comes with a proof of correctness, it doesn't mean it's the correct algorithm for the job.
2. Students will understand that different models / parameters to a model result in different outcomes, and thus even computational problems are not purely computational.
3. Doing theoretical computer science / mathematics is not devoid of ethical considerations.

I would very much appreciate thoughts on the above.

# Answer
> 3 votes
Your understanding of FASTA format is about right. The type of basic problem you're eluding to we term "sequence alignment"- edit distance might be okay for teaching but in practise we use other algorithms, e.g. you might be interested in the Needleman–Wunsch or Smith–Waterman algorithms. Richard Durbin et al. wrote a great book that covers these any much more https://www.amazon.com/Biological-Sequence-Analysis-Probabilistic-Proteins/dp/0521629713

Tracking mutations etc. requires more than just alignment though, see "phylogenetics" (i.e. building genetic trees) and "variant calling". Also checkout what the nextstrain team is doing https://nextstrain.org/ncov/global

In general, looking for practical applications for your algorithms is great, but be very careful before drawing any real world conclusions about the coronavirus outbreak from such analyses

# Answer
> 2 votes
The correct method is laborious and it is better to give the students a ready made tree at GISAID to investigate the dispersion of COVID-19 into Europe.

However, a quick approach to getting to the alignment and drawing a tree is easy and will readily complement your established teaching method. What this will give you is a very different phylogeny from `edit distances` and you would explain the matrix differences between the approaches. I suspect NCBI use Jukes Cantor distances.

1. Go to blast.ncbi.nlm.nih.gov
2. Select nucleotides
3. Paste "MN988669" and hit return
4. This will produce 100 hits
5. On the menu page select the "draw tree option"
6. This will produce a reasonable nucleotide tree for 100 COVID-19 sequences
7. Select "neighbor-joining" rather than "minimum evolution" (this is hierarchical clustering geared to towards mutation rate heterogeneity between taxa)
8. There are various point and click options including "examine alignment"
9. You can then select "minimum evolution" and see the changes in the tree (and there are changes) - TEACHING OBJECTIVE 1 & 2
10. If you want to examine COVID-19 from the European perspective then select a European isolate, e.g. French isolates and input this into the blast. However, European isolates are basically Wuhan original.

I have provided an example below, which I "rerooted" and represents the closest 100 to you sequence using collapsed clade format (your students can undo this, to examine the contents of a "collapsed clade"). The tree shows the dispersion of strains from the Wuhan seafood market.

There's lots of flexibility and the students can do everything inside ½ hour easily and this would complement your approach. The advantage of my approach to teaching phylogeny and regardless of what you do and how you do it Blast is central to rapidly obtaining alignment data for students and investigators alike. We use different blast options but blasting through is prerequisite to understand the diversity and clean some information on the population structure.

---

> Also, are there any recommended data sets where we could do something like track mutations in the virus (and see, e.g., that the NYC outbreak originated from Europe and not China)?

Yes there are, today this would be Mesquite available here. In my opinion this is a bit advanced. You can track amino acid mutations, which in my view is easier.

---
Tags: phylogenetics, sequence-alignment, phylogeny, fasta, covid-19
---


# Question
Title: Transcribing DNA exons and append transcripts into a list
I am trying to find a way to read and transcribe a list of DNA sequences (list of lists) only when the for loop finds a start codon (triplet of the list items) and until it finds a stop codon, over and over, but the if function will only transcribe the first codon and won't go through the whole sequence. I'm trying to do this with a mix of normal python and biopython functions cause I think it's simpler.

This is an example:

Start codon: ATG

Stop codon: TAG

gene\_1= 'ACGGACTATTC'

gene\_2= 'GGCC**ATG**AGTAACGCA**TAG**GGCCC

gene\_3=GGGCCC**ATG**ACGTAC**TAG**GGGCCC**ATG**CATTCA**TAG**

So, the first gene does not contain any transcript, whereas the second contains 1 and the third contains 2. I'm trying to get rid of everything outside these reading frames and append these transcripts into a list that should look something like this.

mRNAs = \['AGUAACGCA', 'ACGUAC', 'CAUUCA'\]

This is what I have so far:

```
genelist=[]

gene_1= 'A','C','G','G','A','C','T','A','T','T','C'
gene_2= 'G','G','C','C','A','T','G','A','G','T','A','A','C','G','C','A','T','A','G','G','G','C','C','C'
gene_3='G','G','G','C','C','C','A','T','G','A','C','G','T','A','C','T','A','G','G','G','G','C','C','C','A','T','G','C','A','T','T','C','A','T','A','G'

genelist.append(gene_1)
genelist.append(gene_2)
genelist.append(gene_3)

def transcription(ORF):
    mRNA= ''
    for i in range(0, len(ORF), 3):
        codon= ORF[i:i+3]
        if codon != 'ATG':
            next(codon)
            if codon == 'ATG':
                mRNA=codon.transcribe()
                if codon == 'TAG':
                    break
    return(mRNA)

mRNAs=[]
for gene in genelist:
    for codon in gene:
        mRNA= transcription(codon)
        mRNAs.append(mRNA)
print(mRNAs)

```

But it is not really giving anything back, I wonder if the code it's too redundant and I really don't need to define a function here, do you know any better way to do this? Thaaanks!!

# Answer
I'm afraid none of your code makes much sense. You have written a function that expects an array, and then you are passing it a single letter

```
for gene in genelist:
    for codon in gene:
        mRNA= transcription(codon)

```

Then, the function also doesn't make sense. What does `next(codon)` mean? How would your script know what a codon is or how to move to the next one? And your `if codon == 'ATG':` is within the `if codon != 'ATG':` block, so it will never be executed by definition.

I think what you are looking for is something like this, using the same basic approach as your script:

```
#!/bin/python3
from Bio.Seq import Seq
from Bio.Seq import transcribe

genelist=[]

gene_1= 'A','C','G','G','A','C','T','A','T','T','C'
gene_2= 'G','G','C','C','A','T','G','A','G','T','A','A','C','G','C','A','T','A','G','G','G','C','C','C'
gene_3='G','G','G','C','C','C','A','T','G','A','C','G','T','A','C','T','A','G','G','G','G','C','C','C','A','T','G','C','A','T','T','C','A','T','A','G'

genelist.append(gene_1)
genelist.append(gene_2)
genelist.append(gene_3)

def transcription(ORF):
    mRNA= ''
    foundStart = False
    foundEnd = False
    for i in range(0, len(ORF), 3):
        codon= "".join(ORF[i:i+3])
        if codon == 'ATG' and not foundStart:
            foundStart = True
        if foundStart and not foundEnd:
            cc=transcribe(codon)
            mRNA = mRNA + transcribe(codon)
        if codon == 'TAG':
            foundEnd = True

    return(mRNA)

mRNAs=[]
for gene in genelist:
    mRNA = transcription(gene)
    mRNAs.append(mRNA)
print(mRNAs)

```

Which returns:

```
$ foo.py
['', '', 'AUGACGUACUAG']

```

Of course, it only returns one mRNA since you only have one gene with an ORF in the 1st frame and your script was only looking at the forward strand and the first frame. In reality, ORFs can be in any of the 6 possible frames. You are also assuming that all ORFs will end with a TAG which is perhaps the case in the species you are working in? I don't remember if there is any genetic code that has TAG as the *only* stop.

If you do actually need to find ORFs in any frame, then I hope the script above should be enough to get you started.
> 0 votes

# Answer
Just in case this *isn't* a homework problem, there are already plenty of programs available to do this. One program I use is `getorf` from the EMBOSS toolkit:

https://emboss.bioinformatics.nl/cgi-bin/emboss/getorf
> 0 votes

---
Tags: python, phylogenetics, homework
---


# Question
Title: Modifying protein_scale in ProteinAnalysis in SeqUtil in Biopython?
Is there a way to give weight arguments to each of the residues when using protein\_scale in the ProteinAnalysis module of Bio.SeqUtils.ProtParam? I thought this would be done with the window and edge arguments, but from the documentation and the question How to allow "half windows" in biopython? I understand that the weights are symmetric and calculated automatically given a window size.

I'd like to understand what the idea behind this was, just to check that I understand this correctly. For example, if one calculates eg. hydrophobicity of a protein, why should it be a value per residue that depends on a particular window?

I was looking for a functionality like where the value given for each residue in the ParamData dictionary would be multiplied by a specified weight, like so:

```Y = ProteinAnalysis('MSTAEFAQLL')
Y.protein_scale(Bio.SeqUtils.ProtParamData.kd,10)
Y.protein_scale(Bio.SeqUtils.ProtParamData.kd,10,edge[0,0.2,0.3,0.1,0.,0.4,0.5,0.8,0.8,1.])```

Of course I can code this myself, but was wondering if that functionality already existed.

# Answer
> 3 votes
There is not a way to control the edge weights as written, but making a modified version is straightforward. Namely make a class that inherits the class `Bio.SeqUtils.ProtParam.ProteinAnalysis` and overwrites the method relevant methods.

Having a gander at the files one sees the following:

```
def protein_scale(self, param_dict, window, edge_list: list):
    """[...]
    """
    # generate the weights
    #   _weight_list returns only one tail. If the list should be
    #   [0.4,0.7,1.0,0.7,0.4] what you actually get from _weights_list
    #   is [0.4,0.7]. The correct calculation is done in the loop.
    weights = self._weight_list(window, edge)
    scores = []
    ...

```

So the weights are generated by the 'private' method `._weight_list`. Underscore means it is for internal use and for production code it is bad to assume that the name of the private method many not change, but Biopython has not changed in years so you'd be safe. And this is for a quick hack anyway...

```
class JankyProteinAnalysis(Bio.SeqUtils.ProtParam.ProteinAnalysis):
    overriding_weights = None

    def _weight_list(self, window, edge):
          if self.overriding_weights is not None: # dodgy override!
               print("Weight method hack done!")
               return overriding_weights
          else: # behave normally
               return super()._weight_list(window, edge)

```

So you'd do:

```
Y = JankyProteinAnalysis('MSTAEFAQLL')
Y.overriding_weights = [0,0.2,0.3,0.1,0.,0.4,0.5,0.8,0.8,1.]
Y.protein_scale(Bio.SeqUtils.ProtParamData.kd,10,float('nan'))

```

Note that the `edge` argument is no longer doing anything.

---
Tags: proteins, python
---


# Question
Title: Bonferroni count = 1 in GO-Slim overrepresentation test (PANTHER) is it possible?
I encountered a problem while performing analysis using the PANTHER website. I run the statistical overrepresentation test for a genelist of 99 genes (GO-Slim BP, Binomial, Bonferroni) in February 2020. My result was Bonferroni count = 1 and a list of several enriched GO terms. Next we run the same analysis with the same settings in July 2020 (two days ago) and the result is dramatically different, the Bonferroni count is 1770 and no enriched terms are found.

I am confused - is the result Bonferroni count = 1 possible? Why these results are so different?

This is how the correction is described on PANTHER website:

> The expression data analysis statistics now include a Bonferroni correction for multiple testing. The Bonferroni correction is important because we are performing many statistical tests (one for each pathway, or each ontology term) at the same time. This correction multiplies the single-test P-value by the number of independent tests to obtain an expected error rate.
>
> For pathways, we now correct the reported P-values by multiplying by the number of associated pathways with two or more genes. Some proteins participate in multiple pathways, so the tests are not completely independent of each other and the Bonferroni correction is conservative. For ontology terms, the simple Bonferroni correction becomes extremely conservative because parent (more general) and child (more specific) terms are not independent at all: any gene or protein associated with a child term is also associated with the parent (and grandparent, etc.) terms as well.
>
> To estimate the number of independent tests for an ontology, we count the number of classes with at least two genes in the reference list that are annotated directly to that class (i.e. not indirectly via an annotation to a more specific subclass).

# Answer
> GO-Slim BP, Binomial, Bonferroni

This seems easy you because I suspect have run too many tests at 5%.

If you run 10 tests at 5% rejection threshold, the Bonferroni correction threshold is 0.5% - you divide the number of tests by the rejection threshold.

One possibility is the rejection threshold is given in number of genems: in one instance you've run a pairwise test with several thousand genes so the rejection threshold has been set accordingly. The other test you ran I suspect was a multi-variate test, thus the shift in rejection threshold measured in number of genes.

---

The altnative is that Bonf = 1 is in fact **P** = 1, in other words there is no significant difference at all. If **P** = 1 then this could be no genes show deviation hence 1770, because 1770/1770 =1

Looking at it again the above seems a better explanation.
> 0 votes

---
Tags: statistics, gene-ontology, go-enrichment
---