# Question
Title: What is the best way to programmatically convert Ensembl ids from one release to other? e.g. transcript ids from release 75 to 100
3 solutions that I had to unfortunately rule out.

**Potential solution #1**

As far as I know, the conventional way is to do this is by using CCDS ids (e.g. using biomart). However, the issue is that, as I see on the CCDS project page, the last release of CCDS ids was made almost 3 years ago (!) and it seems to be compatible with Ensembl release 92. Then how should one convert ids to the latest Ensembl release e.g. 100?

| CCDS Release | Species | NCBI Annotation Release | Ensembl Annotation Release | Assembly Name | Assembly ID | Made Public |
| --- | --- | --- | --- | --- | --- | --- |
| 22 | Homo sapiens | 109 | 92 | GRCh38.p12 | GCF\_000001405.38 | 06/14/2018 |

**Potential solution #2**

I am also aware of Ensembl's ID History Converter. However, I am looking for a programmatic approach. Also, I noticed that if I try to convert ~\>20K ids with ID History Converter, in my experience, there is always some issue. e.g. an error I got when I tried to convert ~70K ensembl transcript ids.

As a side note, I wonder why this ID conversion functionality is not available via Ensembl's REST API. To me, id conversion seems to be a basic enough task that almost every bioinformatician has to do. So it would be immensely helpful to be someday able to do this programmatically using Ensembl's mighty a REST API.

**Potential solution #3**

Finally, as another potential solution, over at Ensembl's FTP website, located within the directory of the release to which I am interested in converting my ids to (release 100, target release), I saw that according to a README file, "ensembl\_archive\_" directory should contain "Data on historical Ensembl releases". So, I wondered if this directory contains a map between ids from different releases. However, without any extra info about column names, contents, file formats etc, I couldn't really figure out how to use those files. So, I had to rule this one out too.

```
Index of /pub/release-100/mysql/ensembl_archive_100/
../
CHECKSUMS                                          26-Mar-2020 14:31                 132
ens_release.txt.gz                                 26-Mar-2020 14:31                1098
ensembl_archive_100.sql.gz                         26-Mar-2020 14:31                 943
release_species.txt.gz                             26-Mar-2020 14:31               39620
species.txt.gz                                     26-Mar-2020 14:31                6669

```

Of note, I saw a similar "mysql" directory over at the directory of the release whose ids I want to convert (source release): `http://ftp.ensembl.org/pub/grch37/release-100/mysql/homo_sapiens_core_100_37/` . However, again, I do not understand its contents without any extra information.

So, finally, I am not sure what is the best way to programmatically convert Ensembl ids from between different releases. As a test case scenario, I want to convert all the protein\_coding transcript ids from GRCh37 (release 75) to GRCh38 (release 100). If anybody has any suggestions about this, please kindly let me know.

# Answer
The ID history converter is available as a Perl script.

This script accesses these MySQL tables, so you could query these tables directly.

Note that the mapping session is only one release at a time, so you would need to run 25 queries to go from release 75 to 100.
> 1 votes

# Answer
Here's one way to do it programmatically using the biomaRt package (your potential solution #1):

```
library(biomaRt)
library(dplyr)

mart37 <- useMart(biomart = "ENSEMBL_MART_ENSEMBL",
                  host    = "grch37.ensembl.org",
                  path    = "/biomart/martservice",
                  dataset = "hsapiens_gene_ensembl")

mart38 <- useMart(biomart = "ENSEMBL_MART_ENSEMBL",
                  path    = "/biomart/martservice",
                  dataset = "hsapiens_gene_ensembl")

#pick all genes on chr22 from build 37 as a test dataset::
out37 <- getBM(attributes = c('hgnc_symbol', 'chromosome_name',
                              'ensembl_gene_id', 'ensembl_transcript_id', 'ccds'),
                  filters = c('chromosome_name'),
                  values = list("22"),
                  mart = mart37) %>% #empty string CCDS ID's
  filter(ccds != "")

#NOTE:: to check which filters can be used to query mart37, run this command:
#View(listFilters(mart37))

out38 <- getBM(attributes = c('hgnc_symbol', 'chromosome_name',
                              'ensembl_gene_id', 'ensembl_transcript_id', 'ccds'),
                        filters = c('ccds'),
                        values = list(out37$ccds),
                        mart = mart38)

```

What this does is query the build 38 DB via the CCDS ID's of build 37. You can include other variables in this `out38` dataframe by simply adding to the attributes argument above. And in order to see all of the potential variables that can be included, you just have to run `View(listFilters(mart37))` or `View(listFilters(mart38))`, depending on which build you are interested in.

Does this solve your problem?
> 0 votes

---
Tags: gene, ensembl, identifiers, conversion
---


# Question
Title: What kind of BLAST do I need to do to accomplish this task?
(I am a complete newbie to bioinformatics, so please bear with me!).

I recently used BLAST to compare a 500 bp nucleotide sequence in *S. cerevisiae* to a bunch of other 500 bp nucleotide sequences in various fungi (these other sequences are homologs of the *S. cerevisiae* gene). Ideally, when I click "BLAST", I want to get nothing back. What I'm trying to double check is the assumption that my *S. cerevisiae* sequence is NOT similar to my other input sequences.

Here's what I ran on the computer, using the BLAST website:

After I run this query, I only get one sequence back:

which tells me that the sequence I got as an output needs to be excluded from my dataset, as it's too similar to the *S. cerevisiae* sequence.

The problem is, I need to run this type of query and analysis about 6000 more times for 6000 different *S. cerevisiae* genes and their homologs. I don't know how to do this programmatically. The Biopython packages seem to only compare one sequence against the whole NCBI database, and I was really struggling to understand how to do this via the command line (was struggling with the package installations a bit, and it also seemed to only be for one sequence input).

Is what I'm doing here a multiple sequence alignment? When I looked into programmatically doing an MSA, the output seemed to only be the alignment, but I want to see an output like the website gives me. I'm not really sure what I should be looking into.

Any guidance would help!

# Answer
## Clarification

Just making sure - you have ~6000 datasets - one for each *S.cerevisiae* gene of interest - which are made up of homologues of each respective gene? And you want to filter out genes that are too similar?

In that case, it would be helpful to know how you're defining a sequence as 'similar enough' for it to be filtered. Given they're all homologues, it's likely there's going to be a high baseline of sequence similarity. In your example, you consider 97.40% identity across the whole query as different enough to be filtered out, so I'm getting the feeling you specifically want to filter out identical sequences?

Also, when you talk about genes here, do you want to filter out gene **products** (i.e. proteins) that are too similar to the actual coding sequences (CDS)? The CDS of two homologues might be X% different, but that difference might not be as high once translated into proteins due to how codons work. Proteins can go down to levels of similarity as low as ~30% identity and still have the same function, so depending on the context of your work (taxonomic, structural, functional) you might want to translate these CDS into protein, so further clarification would be helpful!

## The actual bioinformatics question

To answer your question of programming, if the context of this comparison is always X *S.cerevisiae* gene vs. Y dataset of homologues, a MSA is definitely not the answer (you would do that if you wanted to see how all of the homologues compare to each other, rather than to specifically X gene). Also, using the web-client for NCBI BLAST will take a lot of time, so that's also not recommended for the kind of throughput you need.

Pairwise sequence alignments that you do on your local machine/server between X genes and all of its respective Y dataset sequences is probably the best way to go. I would suggest looking into **needleall** \- it allows you to compare X sequences to Y sequences in a pairwise way, though it performs a **global alignment** rather than **local** like BLAST, so it might be slower.

I don't know what your computational resources are, so in case that might be too computationally heavy, the faster alternative would be to install the **BLAST+ suite**, which allows you to run BLAST queries like you did on the web browser, except on your own machine. Specifically, you could make your own custom individual BLAST databases out of each Y dataset using `makeblastdb` from the suite. Then you use `blastn` (or `blastp` if you decide to translate your sequences) on each dataset with its respective S.cerevisae sequence.

BLAST+ also allows you to customise the actual output as well, making it easier to parse your results programatically for when you want to remove sequences that are too similar. There are many examples online of how to do this, but have a go at it on your own first with the BLAST+ documentation I linked :)
> 6 votes

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


# Question
Title: Nextflow: run sbatch (slurm job script) inside a process
I was trying to run a slurm job script inside nextflow process. Is it possible? For instance, I have a job script 'PRE\_PHASING\_QC.JOB' and outside nextflow it runs perfectly, I want to integrate it in a nextflow script. May be my example code could be an explanation of what I want.

```
prephase='PRE_PHASING_QC.JOB'
process prephasing {
input:
file(prephasing) from prephase
"""
sbatch "${prephasing}
"""
}

```

The error:

```
sbatch: error: This does not look like a batch script
...

```

I tried to make my job script executable, put it in bin and then include it in the nextflow process according to @Steve. So, the script runs without any error but got no output. The original job script takes 3 files (bed,bim,fam) as input and produce multiple outputs. I guess that I have to declare outputs inside nextflow process. So far I tried as follows:

```
process prephasing {
publishDir params.out, mode: 'copy'
output:
tuple val(chr),val(num),val(steps),path("${chr}${num}.${steps}.{bed,bim,fam,log}") into out1
script:
"""
PRE_PHASING_QC.job
"""

```

It gives:

```
Error executing process > 'prephasing'
Caused by:
Missing value declared as output parameter: chr

```

My expected output would be chr{1-22}.step{1-10}.{bed,bim,fam,log}

```
#!/bin/bash
#SBATCH --job-name=prephase
#SBATCH --account=p697_norment
#SBATCH --time=1:00:00
#SBATCH --ntasks=1
#SBATCH --mem-per-cpu=7600M
#SBATCH --cpus-per-task=7
#SBATCH --array=1-22

source /cluster/bin/jobsetup

set -o errexit

module load plink

export prefix=chr${SLURM_ARRAY_TASK_ID}
export hrc_sites_file=/cluster/projects/p697/projects/moba_qc_imputation/resources/HRC/HRC.r1-1.GRCh37.wgs.mac5.sites.tab.gz
export HRC_1000G_check_bim=/cluster/projects/p697/github/norment/moba_qc_imputation/software/HRC-1000G-check-bim.pl

# Step 1a (exclude duplicated), 1b (exclude indels and non-ACTG alleles), 1c (exclude ambiguous SNPs)

plink --bfile hce-ec-eur-fin-pass-qc --chr ${SLURM_ARRAY_TASK_ID} --out $prefix --make-bed || true
plink --bfile $prefix --list-duplicate-vars ids-only suppress first --out $prefix.stepla || true
plink --bfile $prefix --exclude $prefix.stepla.dupvar --make-bed --out $prefix.stepla || true
plink --bfile $prefix.stepla --snps-only just-acgt --make-bed --out $prefix.steplb || true

awk '{ if (($5=="T" && $6=="A")||($5=="A" && $6=="T")ii($5=="C" && $6=="G")||($5=="G" && $6=="C")) print $2, "ambig" ; else print $2 ;}' $prefix. steplb.bim | grep ambig > $prefix.steplb.snplist.txt || true
plink --bfile $prefix.steplb --exclude $prefix.steplb.snplist.txt --make-bed --out $prefix.steplc || true

# Step 2. Zero out Mendelian errors
plink --bfile $prefix.steplc --me 1 1 --set-me-missing --make-bed --out $prefix.step2 || true

# Step 3. Missingness and HWE + Step 4. Compute allele frequencies
# NB! --mind 0.02 filter is removed from the next command, to avoid sample QC per chunk, which may leave inconsistent number of individuals avail able per chunk.
plink --bfile $prefix.step2 --hwe le-6 --geno 0.02 --make-bed --out $prefix.step3 || true
plink --bfile $prefix.step3 --freq --make-bed --out $prefix.step4 || true

# Step 4-10. Run Oxford perl script and apply first 5 lines of its Run-plink.sh command, followed by '--update-name' command (update marker name s)
perl ${HRC-IRC_1000G_check_bim} -b $prefix.step4.bim -f $prefix.step4.frq -r ${hrc_sites_file} -h || true
plink --bfile $prefix.step4 --exclude Exclude-$prefix.step4-HRC.txt --make-bed --out $prefix.step5 || true
plink --bfile $prefix.step5 --update-chr Chromosome-$prefix.step4-HRC.txt --make-bed --out $prefix.step6 || true
plink --bfile $prefix.step6 --update-map Position-$prefix.step4-HRC.txt --make-bed --out $prefix.step7 || true
plink --bfile $prefix.step7 --flip Strand-Flip-$prefix.step4-HRC.txt --make-bed --out $prefix.step8 || true
plink --bfile $prefix.step8 --a2-allele Force-Allele1-$prefix.step4-HRC.txt --make-bed --out $prefix.step9 || true
plink --bfile $prefix.step9 --update-name ID-$prefix.step4-HRC.txt --real-ref-alleles --make-bed --out $prefix.step10 || true

```

# Answer
Using Nextflow to submit jobs this way makes little sense, since:

> sbatch exits immediately after the script is successfully transferred to the Slurm controller and assigned a Slurm job ID. The batch script is not necessarily granted resources immediately, it may sit in the queue of pending jobs for some time before its required resources become available.

Assuming your job script is basically a shell script, make it executable (e.g. `chmod +x PRE_PHASING_QC.JOB`) and move it into a folder called `bin` in the root directory of your project repository. Nextflow automatically adds this folder to your $PATH, making scripts in this folder accessible in your pipeline without needing localize them and without needing to specify an absolute path to invoke them. Your code then becomes:

```
process prephasing {

    """
    PRE_PHASING_QC.JOB [options] <args>
    """
}

```

Finally, add `process.executor = 'slurm'` to your `nextflow.config` file. You may need to configure some other options, like resource requests, but that should get your script integrated into your pipeline without much refactoring, if any.

---

---

Contents of `PRE_PHASING_QC.sh`:

```
#!/bin/bash

set -euo pipefail

hrc_sites_file="${1}"
input_prefix="${2}"
chrom="chr${3}"

# Step 1a (exclude duplicated), 1b (exclude indels and non-ACTG alleles), 1c (exclude ambiguous SNPs)

plink --bfile "${input_prefix}" --chr "${3}" --out "${chrom}" --make-bed
plink --bfile "${chrom}" --list-duplicate-vars ids-only suppress first --out "${chrom}.step1a"
plink --bfile "${chrom}" --exclude "${chrom}.step1a.dupvar" --make-bed --out "${chrom}.step1a"
plink --bfile "${chrom}.step1a" --snps-only just-acgt --make-bed --out "${chrom}.step1b"

awkcmd='
($5 == "T" && $6 == "A") || ($5 == "A" && $6 == "T") ||
($5 == "C" && $6 == "G") || ($5 == "G" && $6 == "C")
{
    print $2, "ambig"
    next
}
{
    print $2
}'

awk "${awkcmd}" "${chrom}.step1b" | grep "ambig" > "${chrom}.step1b.snplist.txt"
plink --bfile "${chrom}.step1b" --exclude "${chrom}.step1b.snplist.txt" --make-bed --out "${chrom}.step1c"

# Step 2. Zero out Mendelian errors
plink --bfile "${chrom}.step1c" --me 1 1 --set-me-missing --make-bed --out "${chrom}.step2"

# Step 3. Missingness and HWE + Step 4. Compute allele frequencies
# NB! --mind 0.02 filter is removed from the next command, to avoid sample QC per chunk, which may leave inconsistent number of individuals avail able per chunk.
plink --bfile "${chrom}.step2" --hwe 1e-6 --geno 0.02 --make-bed --out "${chrom}.step3"
plink --bfile "${chrom}.step3" --freq --make-bed --out "${chrom}.step4"

# Step 4-10. Run Oxford perl script and apply first 5 lines of its Run-plink.sh command, followed by '--update-name' command (update marker name s)
HRC-1000G-check-bim.pl -b "${chrom}.step4.bim" -f "${chrom}.step4.frq" -r "${hrc_sites_file}" -h
plink --bfile "${chrom}.step4" --exclude "Exclude-${chrom}.step4-HRC.txt" --make-bed --out "${chrom}.step5"
plink --bfile "${chrom}.step5" --update-chr "Chromosome-${chrom}.step4-HRC.txt" --make-bed --out "${chrom}.step6"
plink --bfile "${chrom}.step6" --update-map "Position-${chrom}.step4-HRC.txt" --make-bed --out "${chrom}.step7"
plink --bfile "${chrom}.step7" --flip "Strand-Flip-${chrom}.step4-HRC.txt" --make-bed --out "${chrom}.step8"
plink --bfile "${chrom}.step8" --a2-allele "Force-Allele1-${chrom}.step4-HRC.txt" --make-bed --out "${chrom}.step9"
plink --bfile "${chrom}.step9" --update-name "ID-${chrom}.step4-HRC.txt" --real-ref-alleles --make-bed --out "${chrom}.step10"

```

Make `PRE_PHASING_QC.sh` executable (`chmod +x PRE_PHASING_QC.sh`), and replace your job script with the following. This file won't be needed by your Nextflow process, but shows you how the above shell script can be invoked:

```
#!/bin/bash
#SBATCH --job-name=prephase
#SBATCH --account=p697_norment
#SBATCH --time=1:00:00
#SBATCH --ntasks=1
#SBATCH --mem-per-cpu=7600M
#SBATCH --cpus-per-task=7
#SBATCH --array=1-22

set -e

. /cluster/bin/jobsetup

export PATH="/cluster/projects/p697/github/norment/moba_qc_imputation/software:${PATH}"

module load plink

cd /path/to/hce-ec-eur-fin-pass-qc

/path/to/PRE_PHASING_QC.sh \
    /cluster/projects/p697/projects/moba_qc_imputation/resources/HRC/HRC.r1-1.GRCh37.wgs.mac5.sites.tab.gz \
    hce-ec-eur-fin-pass-qc \
    "${SLURM_ARRAY_TASK_ID}"

```

Copy `PRE_PHASING_QC.sh` and the Oxford Perl script, `HRC-1000G-check-bim.pl`, into your project's `bin/` directory. You may also need to make the Perl script executable. Then your Nextflow process becomes:

```
hrc_sites_file = file(params.hrc_sites_file)

Channel
    .fromFilePairs( "${params.input_dir}/*.{bed,bim,fam}", size:3 )
    .set { plink_data }

Channel
    .of(1..22)
    .set { chroms }

process prephasing {

    input:
    tuple val(prefix), path(bfiles) from plink_data
    each chrom from chroms
    path hrc_sites_file

    """
    PRE_PHASING_QC.sh "${hrc_sites_file}" "${prefix}" "${chrom}"
    """
}

```

Once the above runs successfully, you can declare your process outputs. Please note that none of the above was tested; hopefully not too many typos, but YMMV.
> 2 votes

---
Tags: nextflow
---


title: Do proteins move inside the cell?

text:
# Question
Title: Do proteins move inside the cell?
Do proteins move inside the cell and how fast do they move or move? Please support your answer with articles or books.

# Answer
> 2 votes
Proteins move in cells or get transported in the cell. You can go to NCBI (https://www.ncbi.nlm.nih.gov/), and start searching for protein dynamics within cells. Also there is a technique called FRAP (https://en.wikipedia.org/wiki/Fluorescence\_recovery\_after\_photobleaching) which is used to measure kinetics and dynamics.
Also, this differs probably a lot between the cellular system you are asking for. You can imagine that some cells have a more active transport system because stretched out (neurons) vs others with maybe less active transport systems, because of a more static nature.

---
Tags: proteins
---


# Question
Title: Toptable error, wont recognize condition
I am getting a rather strange error from topTable. When I run my code I get `Error in fit$coefficients[, coef] : subscript out of bounds` from topTable as if it is not recognizing the condition I am asking it to return the top ten DEGs of. I have provided a bit more context and some descriptions of my variables below. Let me know if you need any more information or if anything is unclear.

**Variables:**

`dge` is a DGEList object of RNA seq data.

`dge$counts` is the raw counts.

`dge$samples$dge_meta` is the metadata and log2 transformed counts.

```
> colnames(dge$samples$dge_meta)
 [1] "Lab.sample.ID"
 [2] "Accession.Number"
 [3] "Alias.used.in.some.figures"
 [4] "Developmental.stage"
 [5] "Infected"
 [6] "Batch"
 [7] "Trimmed"
 [8] "Number.of.reads.that.pass.Illumina.filter"
 [9] "Total.number.of.reads.mapped"
[10] "%.of.total.reads.mapped"
[11] "Reads.mapped.to.T..cruzi.Esmeraldo.haplotype"
[12] "Reads.mapped.to.hg19"
[13] "%.of.mapped.reads.belonging.to.T..cruzi.Esmeraldo.haplotype"
[14] "%.of.mapped.reads.belonging.to.hg19"
[15] "sample"
[16] "hpi"

```

Lastly `keep_it` is a binary vector used to filter our zero count genes.

**Code:**

```
mm <- model.matrix(~~hpi*Infected, data = dge$samples$dge_meta) #making our design matrix

more_than_1cpm_model <- voom(dge$counts[keep_it,], mm, plot = TRUE) #this one is filtered by the 25% greater than 1 cpm

cutoff <-  1e-05

empBayesfit <- eBayes(fit_gt1cpm)

topTable(empBayesfit,
         coef = c("Infected"),
         p.value = cutoff, number = 10)

```

**Error message:**

```
Error in fit$coefficients[, coef] : subscript out of bounds

```

To understand the context better I ran rlang::last\_error and rlang::last\_trace as well.

```
> rlang::last_error()
<error/dplyr_error>
arrange() failed at implicit mutate() step.
* Problem with `mutate()` input `..1`.
x Input `..1` must be a vector, not a `standardGeneric` object.
i Input `..1` is `samples`.
Backtrace:
Run `rlang::last_trace()` to see the full context.
> rlang::last_trace()
<error/dplyr_error>
arrange() failed at implicit mutate() step.
* Problem with `mutate()` input `..1`.
x Input `..1` must be a vector, not a `standardGeneric` object.
i Input `..1` is `samples`.
Backtrace:
     x
  1. +-dplyr::arrange(x, samples)
  2. +-dplyr:::arrange.data.frame(x, samples)
  3. | \-dplyr:::arrange_rows(.data, dots)
  4. |   +-base::withCallingHandlers(...)
  5. |   +-dplyr::transmute(new_data_frame(.data), !!!quosures)
  6. |   \-dplyr:::transmute.data.frame(new_data_frame(.data), !!!quosures)
  7. |     +-dplyr::mutate(.data, ..., .keep = "none")
  8. |     \-dplyr:::mutate.data.frame(.data, ..., .keep = "none")
  9. |       \-dplyr:::mutate_cols(.data, ...)
 10. |         +-base::withCallingHandlers(...)
 11. |         \-mask$eval_all_mutate(quo)
 12. +-dplyr:::abort_glue(...)
 13. | +-rlang::exec(abort, class = class, !!!data)
 14. | \-(function (message = NULL, class = NULL, ..., trace = NULL, parent = NULL, ...
 15. |   \-rlang:::signal_abort(cnd)
 16. |     \-base::signalCondition(cnd)
 17. +-(function (e) ...
 18. | \-rlang::abort(...)
 19. |   \-rlang:::signal_abort(cnd)
 20. |     \-base::signalCondition(cnd)
 21. \-(function (cnd) ...

```

# Answer
The error is telling you that there isn't any Infected column on the model.matrix. Check the column names of `mm`. Also not sure that having two `~` will work well.
> 1 votes

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


# Question
Title: Annotating gene names or gene IDs to a dataframe containing SNPs?
I have a large data-frame (excel file) of SNPs with genotyping data. I need to filter the data to get SNP information of a specific gene alone. The list is too long to do it manually. I was wondering of annotating the SNPs with the gene names so that I can filter using the gene name. Does anyone know of a good R package to perform the task?

# Answer
I assume your SNP file contains genomic positions - please clarify if not. For example, you could have an rsID column, chromsome, and position column that contain "rs3088379", "chr1", "54225281".

**Approach #1:** use R, as you suggest.

`ensembldb` is an R package that may be what you're looking for: https://bioconductor.org/packages/release/bioc/html/ensembldb.html.

You could map the genomic positions of your SNPs to gene names by adapting their workflow: https://bioconductor.org/packages/release/bioc/vignettes/ensembldb/inst/doc/coordinate-mapping.html. After you've mapped the genomic position, you could filter for your gene name.

**Approach #2:** use bedtools.

You can convert your Excel file of genomic positions to a tab-delimited BED format. https://bedtools.readthedocs.io/en/latest/content/overview.html has an overview. So, your first SNP would turn into

```
chr1 54225280 54225281 rs3088379

```

--\> make sure you turn 1-based to 0-based coordinates! You could then download the coordinates of your gene in BED format from your favorite source, like UCSC Table Browser.

Then, you run `bedtools intersect` with your SNP BED and your gene BED. What remains is just that gene's SNPs.
> 2 votes

---
Tags: r, gene, sequence-annotation, snp, filtering
---


# Question
Title: File format of substitution matrix in clustalw
I need to set the substitution matrix used by command line CLUSTALW when comparing DNA sequences to:

```
 0 -1 -1 -1
-1  0 -1 -1
-1 -1  0 -1
-1 -1 -1  0

```

from my understanding I need to pass the location of a file with the matrix to the `-dnamatrix` parameter, however, I cannot find out what format the file is supposed to be in. Any idea?

# Answer
This example seems to have information on this, looks like it should have this format:

```
   A   G   C   T   U   *
A  1   -1  -1  -1  -1  -1
G  -1  1   -1  -1  -1  -1
C  -1  -1  1   -1  -1  -1
T  -1  -1  -1  1   -1  -1
U  -1  -1  -1  -1  1   -1
*  -1  -1  -1  -1  -1  -1

```

However, the page I linked is specifically about a bug regarding ClustalW **ignoring** custom made matrices. It's a pretty old post so it might have been fixed since. However, if you encounter that bug, that same page thankfully does contain a workaround!
> 3 votes

---
Tags: sequence-alignment, multiple-sequence-alignment, substitution-model
---


# Question
Title: Understanding ViennaRNA RNAdistance scoring table
I'm trying to compare the output of 2 different algorithms of RNA structure prediction (my implementation of Nussinov vs RNA-mfold algorithm) using the RNAdistance algorithm that is part of ViennaRNA package. I'm getting a distance score of 38 between Nussinov prediction and the actual NMR structure, and a score of 14 between the m-fold algorithm and the NMR structure and I'm trying to understand whether it is a big difference.

In order to do that I have to look at the scoring table that is used by RNAdistance and found the following `score table` in the source code:

```
RIVATE CostMatrix  UsualCost =
{

/*    Null,       U,        P,        H,        B,        I,        M,        S,        E,        R     */

   {        0,        1,        2,        2,        2,        2,        2,        1,        1, DIST_INF},   /* Null replaced */
   {        1,        0,        1, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF},   /* U    replaced */
   {        2,        1,        0, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF},   /* P    replaced */
   {        2, DIST_INF, DIST_INF,        0,        2,        2,        2, DIST_INF, DIST_INF, DIST_INF},   /* H    replaced */
   {        2, DIST_INF, DIST_INF,        2,        0,        1,        2, DIST_INF, DIST_INF, DIST_INF},   /* B    replaced */
   {        2, DIST_INF, DIST_INF,        2,        1,        0,        2, DIST_INF, DIST_INF, DIST_INF},   /* I    replaced */
   {        2, DIST_INF, DIST_INF,        2,        2,        2,        0, DIST_INF, DIST_INF, DIST_INF},   /* M    replaced */
   {        1, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF,        0, DIST_INF, DIST_INF},   /* S    replaced */
   {        1, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF,        0, DIST_INF},   /* E    replaced */
   { DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF,        0},   /* R    replaced */

};

```

I also found the enter link description herefollowing explanation, but **I do not understand what is the meaning of the table letters: U,P,H,B,I,MS,E,R ?**

What is the meaning of the score of \[H,B\] = 2 for example? If anyone can explain the table it will be a great help in my project, thank you!

# Answer
> 2 votes
Ok, the answer was right here and I missed it, posting it in case anyone else will do the same mistake as me in the future.

```
U = unpaired
P = paired
H = hairpin
B = bulge
I = internal loop
M = multiloop
S = stem (or stack)
E = external loop

```

---
Tags: sequence-alignment, rna, rna-alignment, rna-structure
---


# Question
Title: How low does an E-Value need to be for a gene to code for the same protein as query?
I'm looking at a bacterial genome and trying to find specific PEGs (Protein Encoding Genes) by BLAST searching. I'm taking the genes found with the lowest E-value in order to psiBLAST them to check it codes for the same protein. I'm not exactly sure of how high the E-Value needs to be for me to conclude that the protein is not coded for by the subject gene.

Does anyone have any advice on how to make this process a little more efficient/scientific? Thank you. (By the way, I'm using BLAST with a BLOUSUM62 matrix.)

# Answer
> 3 votes
I assume that you mean "is homologous to" by "encodes the same protein as". Identical ("same") proteins should be fairly trivial to identify using a tool such as tblastx.

You should probably not use E-value for this. E-value is similar to a p-value, in that it represents some aspect of the data (in this case a sequence and a database) compared to some null hypothesis (meaning in this case statistics on the database). Therefore, it depends both on the specifics of the data and the hypothesis, and there is no universal threshold that will work across all sequences and all databases. (See discussion of the meaning of e-value.)

E-value does not in any way represent the content of the sequence or its relationship to a sequence in the database. You could have E-values close to zero or E-values \> 10 for a perfect sequence match or a very loose sequence match, depending on the input sequence and the database.

However, BLAST output does include other information that might be helpful. You could (for example) look at the length of the alignment and the number of mismatches/identities in the alignment. This will give you a pretty good guess at how similar the two sequences are; unless they are very divergent you would expect them to align across most of the length with somewhat high protein identity or similarity (\>25% probably?). This process will still not work 100% of the time, but I think it is closer to your intent than any given approach involving E-values.

---
Tags: blast
---


# Question
Title: Why did expression based subtypng of breast cancer gain much more acceptance than others
This is may not be entirely technical question but rather a academic question. But the technique behind the application is within the scope of bioinformatics. So I would try to ask here that:

In each cancer type, there have been tons of papers that describe RNA-Seq based subtyping resulting in molecularly distinct subgroups with survival difference. However, only very few of them gained acceptance at clinical setting, such as breast cancer (luminal subtype A, luminal subtype B, normal breast-like subtype, HER-2 overexpression subtype and basal-like subtype).

My question is: why did this molecular subtyping of breast cancer gain much acceptance but not others? Did this one own some advantages that other did not?

# Answer
> 6 votes
## In two words: incidence and funding

I'm not an expert on this topic, but I assume it has something to do with the incidence of breast cancer itself:

> Breast cancer is the most common cancer in American women, except for skin cancers. Currently, the average risk of a woman in the United States developing breast cancer sometime in her life is about 13%. This means there is a 1 in 8 chance she will develop breast cancer. This also means there is a 7 in 8 chance she will never have the disease.

This incidence most likely amplifies how big of a problem society/funders/governments consider breast cancer to be. For example, a study performed by Northwestern University’s Feinberg School of Medicine and Robert H. Lurie Comprehensive Cancer Center showed that breast cancer research is by and large one of the most well-funded areas of cancer research.

This higher amount of funding leads to more and better resources and research into every aspect of breast cancer, from identification to treatment. I would wager better funding has had the effect of increasing the level of research performed in the usage of sub-typing for identifying people at risk of getting breast cancer. **Given positive results and lots of literature backing up such a method as useful, its levels of acceptance in clinical settings increased.**

---
Tags: rna-seq, gene-expression
---


# Question
Title: How to obtain NCBI Nucleotide (Entrez) IDs for all E coli sequences?
I would like to gather a list of ALL E coli NCBI (Entrez/Nucleotide) IDs so that I can download fna and faa files, through a function like this:

```
# python code #
from Bio import Entrez
from glob import glob

def dl_genome(id, folder='genomes'):
    files=glob('%s/*.fna'%folder)
    out_file = '%s/%s.fna'%(folder, id)

    if out_file in files:
        print (out_file, 'already downloaded')
        return
    else:
        print ('downloading %s from NCBI'%id)

    Entrez.email = ""
    handle = Entrez.efetch(db="nucleotide", id=id, rettype="fna", retmode="text")
    with open(out_file, 'w') as fout:
        fout.write(handle.read())

for strain in targetStrainIDs:
    dl_genome(strain, folder='genomes')

```

My problem is that I don't have a good way to get the E coli strain names and their corresponding NCBI IDs for Entrez. Is there a way to get these NCBI IDs without manually searching for them?

**Context**

I have been trying to develop an E coli pan-genome for one of my projects. I initially just downloaded the relevant genomes (fna, faa, & features.tab files) from the PATRIC database. However, I've come to realize that the reference sequence on PATRIC has issues. Specifically, the reference sequence on PATRIC's fna/faa files for E coli K12 MG1655 notably do NOT list genes like b2031 or b2033, despite both being present in this strain.

The NCBI reference genome DOES include these genes (with proper annotations), so I'd like to just download what I need from NCBI instead. However, I'd like a systematic way to download all the fna/faa files since my pan-genome is eventually going to be built from hundreds of of these fna/faa files.

Any help would be appreciated.

# Answer
You should check out NCBI Datasets for this sort of thing. Since this is bacteria, you won't be able to use the web interface but there's an API and a command-line tool that you can use.

For example, you can use the `datasets` command line tool to download the E. coli RefSeq genome assemblies as follows:

```
## see how many genomes are in scope
$ datasets summary genome taxon 562 --limit none --refseq
{"total_count":22287}

## download data excluding gff3, rna and genome seq files
## as an example but drop --exclude-seq to download genome seq as well
$ datasets download genome taxon 562 --refseq --exclude-gff3 --exclude-rna --exclude-seq --filename ecoli_genomes.zip

```
> 1 votes

# Answer
Do you mean the NCBI taxonomy IDs? If so, I've personally always preferred the way UniProt displays NCBI taxonomy, and it's easier to download IDs in mass.

For example, if you find *E.coli*'s entry: https://www.uniprot.org/taxonomy/562

You can then click on "all lower taxonomy nodes", which will show all *E.coli* strains. If you click 'Preview first 10' before downloading, you'll see that the first column is the NCBI taxonomy ID of each strain of *E.coli*.

For example, *E.coli* strain 0.1197 is reported to have the ID **1005443** from this preview page. Indeed, if we search this ID on NCBI's taxonomy page, we do get the correct entry: https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=1005443

If this is what you want, click 'download all' to get a nice tab-separated file of all the strains of *E.coli* and their NCBI taxonomy IDs. This file could then be read and processed before you send API requests to Entrez.

Hope this helps!
> 6 votes

---
Tags: phylogenetics, entrez
---


# Question
Title: Differential gene expression analysis of time series with replicates
I have a dataset that has two groups, perturbed vs control. Each group has 3 replicates. Each replicate has 8 time points. How do we do Differential gene expression analysis to find significantly expressed genes? Any packages in R or web tools etc?

# Answer
R packages like DESeq2, edgeR, and sleuth will help you perform differential gene expression analysis. These are typically used for control-vs-treatment experimental designs. So you could use them like this and perform a bunch of pair-wise comparisons between the time points (i.e. what are the differentially expressed genes between timepoint $t\_0$ and timepoint $t\_1$, $t\_0$ and $t\_2$, etc.).

But all of these methods make use of generalized linear models to perform the statistical inference on coefficients. For the simple control-vs-treatment design, above, this looks something like

$$g \left(\mathbb{E} \left\[Y\_i | X \right\] \right) = \mu\_i + \beta\_i X$$

where $g$ is the link function, $\mu\_i$ is the mean expression of gene $i$, $\beta\_i$ is the expression fold change estimate, and $X$ is the indicator variable saying if the sample is a treatment or control sample.

You can modify this equation by adding a time component to it.

$$g \left(\mathbb{E} \left\[Y\_i | X \right\] \right) = \mu\_i + \beta\_i^{(group)} X + \beta\_i^{(time)} t$$

If you specify the timepoint of the sample in the experimental design matrix for DESeq2, edgeR, or Sleuth, you can get information about the effect of the perturbation as well as the effect over time.

You can try more complicated designs, too, like a coefficient for the interaction between timepoint and perturbation status, but with 3 replicates in each, it's probably best to stick with simple, linear models.
> 1 votes

---
Tags: rna-seq, gene, differential-expression
---


# Question
Title: How to get 2d/3d structural coordinates from a PDB file?
Is there a way we can get 2d/3d structural coordinates like in (x/y/z) manner from a pdb file.

I have looked a believe that the coordinates are in below line:

```
ATOM   1067  NH1 ARG A 141     -10.147   7.455  -6.079  1.00 23.24           N

```

but i am unable to make a sense from it as i don't know which one is (x/y/z) in (2d/3d) coordinate system.

If anyone describe it, that'll be a great help.

Thanks

# Answer
> 3 votes
Molecular structures are 3D. The records (lines) in the PDB files are in so-called fixed column format.
The fields don't need to be separated with blank characters, but they must be in specified columns. For example, the X coordinate is in columns 31-38. Here is the specification of the ATOM record from PDB v3.3.

```
COLUMNS        DATA  TYPE    FIELD        DEFINITION
-------------------------------------------------------------------------------------
 1 -  6        Record name   "ATOM  "
 7 - 11        Integer       serial       Atom  serial number.
13 - 16        Atom          name         Atom name.
17             Character     altLoc       Alternate location indicator.
18 - 20        Residue name  resName      Residue name.
22             Character     chainID      Chain identifier.
23 - 26        Integer       resSeq       Residue sequence number.
27             AChar         iCode        Code for insertion of residues.
31 - 38        Real(8.3)     x            Orthogonal coordinates for X in Angstroms.
39 - 46        Real(8.3)     y            Orthogonal coordinates for Y in Angstroms.
47 - 54        Real(8.3)     z            Orthogonal coordinates for Z in Angstroms.
55 - 60        Real(6.2)     occupancy    Occupancy.
61 - 66        Real(6.2)     tempFactor   Temperature  factor.
77 - 78        LString(2)    element      Element symbol, right-justified.
79 - 80        LString(2)    charge       Charge  on the atom.

```

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


# Question
Title: Input from multiple channel, merge and use in nextflow
I need to merge multiple channel to reuse them in two processes. So far I tried as follows:

```
channel.from([15,20,21]).into {k1;k2}
channel.from(['a','b','c']).into {m1;m2}
channel.from([1,2,3]).into {l1;l2}
para1=m1.merge(k1).merge(l1)
para2=m2.merge(k2).merge(l2)
process p1{
input:
each m, k, l from para1
"""
echo $m, $k, $l
"""
}
process p2{
input:
each m, k, l from para2
"""
echo $m, $k, $l
"""
}

```

The error I am getting:

```
Unknown process directive: '_in_each'

Did you mean of these?
    each

```

If I use para1.view() I can see the merged elements as expected. Any help

# Answer
> 2 votes
I don't think you can unpack tuples using the each qualifier like that. You can, however, unpack them in the script block:

```
process test {

    input:
    each triple from para1

    script:
    def (m, k, l) = triple

    """
    echo "${m} ${k} ${l}"
    """
}

```

Consider also that the merge operator is deprecated and will be removed from future Nextflow releases. Instead, you may need to join by index, for example:

```
def indexedChannel( items ) {
    return Channel.from( items.withIndex() ).map { item, idx -> tuple( idx, item ) }
}

k_ch = indexedChannel( [ 15, 20, 21 ] )
m_ch = indexedChannel( [ 'a', 'b', 'c' ] )
l_ch = indexedChannel( [ 1, 2, 3 ] )

m_ch
    .join( k_ch )
    .join( l_ch )
    .map { idx, m, k, l -> tuple( m, k, l ) }
    .into { para1; para2 }

```

---
Tags: nextflow
---


# Question
Title: How to represent random sequence elements in SBOL?
In our lab, we work with synthetic biology components with partially random sequences (similar to work in directed evolution). So, for example, we have a plasmid design with several components, including one protein coding sequence that will be created from a wild-type sequence via mutagenic PCR.

At the design stage, we don't know the exact DNA sequence for any of the resulting variants of that protein coding sequence. We just know the wild-type sequence and the average number of mutations we'd like to apply. How can we represent that design in SBOL (via the `Component` class)?

Then, when we implement the design, we end up with a single sample that contains plasmids with 10^5-10^6 (or more) different protein coding sequences. We then measure the complete DNA sequence for every plasmid variant, but I don't want to end up with \>10^5 different SBOL definitions to describe the sample. How can we best represent that sample in SBOL (via the `Implementation` class)?

Finally, we also end up measuring the performance of every plasmid variant, and we identify the specific protein coding sequences for useful and/or interesting variants. We then design and build specific, clonal plasmids with those useful/interesting protein coding sequences. Representing those clonal plasmids as SBOL `Components` seems pretty straight-forward, but is there a good way to capture the relationships with the original random-library plasmid design and the specific implementation of that design that we built and measured?

# Answer
> 3 votes
This sounds like a job for the `CombinatorialDerivation` class, which is designed for representing libraries and the relationship between derived designs and the originals.

### About `CombinatorialDerivation`

In SBOL 3, a `CombinatorialDerivation` takes a `Component` in its `template` property to represent the system that you intend to vary:

* Some of the `Features` of the template `Component` are designated as the variables of the library via `VariableFeature` child objects of the `CombinatorialDerivation`, which point at the variables with their `variable` property.
* The values that can be filled in for each variable are then designated via the `variant*` properties of the `VariableFeature`.
* Finally, the distribution of values is determined by the `strategy` property of the `CombinatorialDerivation` and the `cardinality` properties of each variable.

### Representing Random Mutagenesis

In the case of random mutagenesis, if I'm understanding your design correctly, you would represent it with a `CombinatorialDerivation` as follows:

* The `template` is `Component` of type DNA with the original sequence. Each site designated for mutagenesis can be designated with a `SequenceFeature`.
* For each of the mutagenesis sites, create a `VariableFeature` in the `CombinatorialDerivation` with the corresponding `SequenceFeature` as its `variable` and a set of `variant` values for each of the mutant options set by your primers. The cardinality is probably `one` (if every site has to be hit by *something* in your protocol) or `zeroOrOne` (if sites might or might not be hit).
* The `strategy` for the `CombinatorialDerivation` is `sample`, since only a few of the many possibilities will actually be realized.

### Representing Samples and Characterized Constructs

Now, when you take it to the lab, you can represent your aliquot with an `Implementation` whose `wasDerivedFrom` points to the `CombinatorialDerivation`, representing the fact that you've got a randomly generated library whose contents are not precisely known.

Finally, once you've got your isolated products represented as `Component`s, there is a recommended best practice to use `wasDerivedFrom` links to connect the derived `Component` and its varied `Feature`s back to the `CombinatorialDerivation` and its template. If you make a `Collection` of your results, the `Collection` should also link to the `CombinatorialDerivation`. This should then capture the relationship between all of the products, their associated data, and the original mutagenesis design.

---
Tags: sbol, variants, library
---


# Question
Title: How to share a collection of sbol designs in SynBioHub
We currently work on a project that will generate a collection of recombinant plasmids targeting over 3500 genes of the cyanobacterium Synechocystis sp. PCC 6803. We would like to share these designs in SynBioHub but we are not sure how best to go about it. The plasmids have been designed and now they are starting to be produced in batches, checked and finally transformed into Synechocystis to induce a single gene deletion by homologous recombination. During the project, we will have:

* Stage 1 – the initial plasmids designs, waiting for production. They are worth sharing as the community could influence the production queue.
* Stage 2 – the assembled plasmids. In some cases they will have different structure than initial designs from stage 1. If nothing else their description should reflect the status change to 'produced'.
* Stage 3 – plasmids are verified by sequencing (also to check the barcodes are correct) and sequences attached.
* Stage 4 – plasmids will be transformed into Synechocystis and checked by PCR to see if i) the target gene has been knocked out and ii) if the strain is segregated (i.e. all copies of the genome in the strain contain the gene deletion). This information should also be recorded with the description.

Our initial plan was:

* Create a collection with initial designs (component definitions) Stage 1
* As stage 2 progresses add new versions of udpated designs
* Stage 3, 4 update plasmids descriptions (synbio desc or notes) and attach additional files (e.g. from sequencing).

However, our current understanding is the public collections are frozen so stage 2,3,4 cannot be accommodated.

We are seeking advice how to plan for sharing this resource and some clarifications.

1. Is there a way of making collection readable by others without freezing its content?
2. What operations are possible on public component definitions, can the free text descriptions/notes be updated or file attached?
3. Can new component definitions or their versions be added to a published collection
4. For 2) 3) who can perform these operations, the collection’s creator or SynBio curators?

Looking forward to input from SynBioHub community

# Answer
There is a balancing act between allowing all edits, which defeats the purpose of having URIs point to static objects, and allowing no edits, which defeats the purpose of providing some updated information about an object. SynBioHub resolves this issue by allowing a limited amount of edits after publication. To answer your questions specifically:

1. Yes, you can share objects by generating a share link OR adding another owner to an object. Sharing can be done by the original owner of the object.
2. There are some mutable fields that any owner of the object can change, these include the description, notes, source, and references in the Details section. You can also add attachments.
3. Yes, when you publish (Make Public) a new object or collection, you have the option of making it public into an existing collection OR a new collection. A person with curator privileges on the SynBioHub instance you are working on can perform this operation.
4. Answered above.
> 9 votes

---
Tags: sbol
---


# Question
Title: Merge multiple channels each emits multiple files in nextflow
I need to merge multiple channels each emits multiple file. So far I did something as below:

```
process prephasing {
publishDir params.out, mode:'copy'
input:
tuple val(prefix), path(bfiles) from plink_data
each chrom from chroms
path hrc_sites_file
output:
tuple val("${chrom}"),path("chr${chrom}*.{bed,bim,fam}") into qc1_chr
"""
PRE_PASHING_QC2.sh "${hrc_sites_file}" "${prefix}" "${chrom}"
"""
}
channel.fromPath('/some/path/genetic_map_chr*_combined_b37.txt').into {shape_map;i
mp_map}
channel.fromPath('/some/path/*hrc.r1-1.ega.grch*.hap.gz').into {shape_ref;
imp_ref}
channel.fromPath('/some/path/*hrc.r1-1.ega.grch37.chr*.legend.gz').into {s
hape_legend;imp_legend}
channel.fromPath('/some/path/*hrc.r1-1.ega.grch37.chr*.samples').into {sha
pe_sample;imp_sample}
qc1_chr.filter{name, bfiles -> bfiles.collect {it.getBaseName()}.every {it.endsWith("step10")}}.into {shape_in;imp_in}

shape_para=shape_in.merge(shape_map).merge(shape_ref).merge(shape_legend).merge(shape_sample)

process phasing {
publishDir params.out, mode:'copy'
input:
each file(penta) from shape_para
output:
path("*step10*") into phased
script:
def (in, map, ref, legend, sample) = penta
"""
shapeit2 -T 32 --input-bed "${in}" --input-map "${map}" --input-ref "${ref}" "${legend}" "${sample}" --duohmm -O "${in}.phased" --seed 54321
"""
}

```

The "prephasing" process works fine and produce expected output. But the "phasing" process does not produce any output nor any error. Any help?

I have tried according to `@Steve` as below:

```
channel.fromFilePairs('/cluster/projects/p697/projects/moba_qc_imputation/resources/1000Genomes/genetic_map_chr*_combined_b37.txt',size:1).map {group_key, file_list -> tuple(group_key.replaceFirst(/^genetic_map_/,""),file_list.first())}.set {genetic_map}

channel.fromFilePairs('/cluster/projects/p697/projects/moba_qc_imputation/resources/HRC/haplegendsample/*hrc.r1-1.ega.grch37.chr*.haplotypes.{hap.gz,legend.gz,samples}',size:3).map{group_key,file_list -> tuple(group_key.replaceFirst(/^*HRC\.r1-1\.ega\.grch37\./,""),file_list)}.set {ref_panel}

qc1_chr.filter{chrom, bfiles -> bfiles.collect {it.getBaseName()}.every {it.endsWith("step10")}}.join(ref_panel).join(genetic_map).into {shape_in;imp_in}

process phasing {
publishDir params.out, mode:'copy'
input:
tuple val(chrom), path(input_bed), path(intput_ref), path(input_map) from shape_in
output:
tuple val("${chrom}"),path("${chrom}.{phased.haps,phased.sample,log}") into phased
script:
def (bed, bim, fam) = input_bed
def (haplotype, legend, sample) = input_ref
"""
shapeit2 --thread "${task.cpus}" --input-bed "${bed}" "${bim}" "${fam}" --input-map "${input_map}" --input-ref "${haplotype}" "${legend}" "${sample}" --duohmm --output-max "${chrom}.phased" --seed 54321
"""
}

```

But getting no output so far, (attached image)

I have tried to check individual channel separately using channel.view(), identified the problems, modified the code and checked again. Using the code below I can see individual channel emits 22 lines \[chr1-22,\[file\_list\]\]. After merging I can see \[chr1-22,\[joined\_file\_list\].

```
channel
    .fromFilePairs( "output3/*.{bed,bim,fam}", size:3 )
    .set { plink_data }
plink_data.filter{chrom, bfiles -> bfiles.collect {it.getBaseName()}.every {it.endsWith("step10")}}.map{group_key,file_list -> tuple (group_key.
replaceFirst(/^*\.step10/,""), file_list)}.set {in_pl}

channel.fromFilePairs('/cluster/projects/p697/projects/moba_qc_imputation/resources/1000Genomes/genetic_map_chr*_combined_b37.txt',size:1).map {
group_key, file_list -> tuple(group_key.replaceFirst(/^*genetic_map_/,""),file_list.first())}.set {genetic_map}
//genetic_map.view()

channel.fromFilePairs('/cluster/projects/p697/projects/moba_qc_imputation/resources/HRC/haplegendsample/*hrc.r1-1.ega.grch37.chr*.haplotypes.{ha
p.gz,legend.gz,samples}',size:3).map {group_key,file_list -> tuple(group_key.replaceFirst(/^.*.hrc\.r1-1\.ega\.grch37\./,""),file_list)}.set {re
f_panel}
//ref_panel.view()

in_pl.join(ref_panel).join(genetic_map).into {shape_in;imp_in}
process phasing {
//publishDir params.out, mode:'copy'
input:
tuple val(chrom), path(input_bed), path(intput_ref), path(input_map) from shape_in
output:
tuple val("${chrom}"),path("${chrom}.{phased.haps,phased.sample,log}") in
to phased
script:
def (bed, bim, fam) = input_bed
def (haplotype, legend, sample) = input_ref
"""
//shapeit="/cluster/projects/p697/projects/moba_qc_imputation/software/shap
eit"
shapeit --thread "${task.cpus}" --input-bed "${bed}" "${bim}" "${fam}" --input-map "${input_map}" --input-ref "${haplotype}" "${legend}" "${samp
le}" --duohmm --output-max "${chrom}.phased" --seed 54321
"""
}

```

But now I get error:

```
No such variable: chrom

```

But I decalred "chrom" in the input

The input bed|bim|fam files from channel from FilePairs works but from a output channel of previous process does not. If I run `phasing` solo and all the inputs from channel from FilePairs it works. But if I run `prephasing` and `phasing` together, the `prephasing` process works but for `phasing` process not output and no error; Not sure what I am doing wrong.

```
hrc_sites_file = file(params.hrc_sites_file)

Channel.fromFilePairs( "${params.input_dir}/*.{bed,bim,fam}", size:3 ).set { plink_data }

Channel.of(1..22).set { chroms }
params.out='output/'

process prephasing {
publishDir params.out, mode:'copy'
input:
tuple val(prefix), path(bfiles) from plink_data
each chrom from chroms
path hrc_sites_file
output:
tuple val("${chrom}"),path("chr${chrom}*.{bed,bim,fam}") into qc1_chr
tuple val("${chrom}"),path("chr${chrom}*log") into qc1_log
tuple val("${chrom}"),path("*-chr${chrom}.step4-HRC.txt") into qc1_hrc
path("chr${chrom}.step1b.snplist.txt") into qc1_snps
path("chr${chrom}.step1a.dupvar") into qc1_dups
"""
PRE_PASHING_QC2.sh "${hrc_sites_file}" "${prefix}" "${chrom}"
"""
}
channel.fromFilePairs('/cluster/projects/p697/projects/moba_qc_imputation/resources/1000Genomes/genetic_map_chr*_combined_b37.txt',size:1).map {group_key, file_list -> tuple(group_key.replaceFirst(/^*gen
etic_map_/,""),file_list.first())}.set {genetic_map}

channel.fromFilePairs('/cluster/projects/p697/projects/moba_qc_imputation/resources/HRC/haplegendsample/*hrc.r1-1.ega.grch37.chr*.haplotypes.{hap.gz,legend.gz,samples}',size:3).map{group_key,file_list ->
 tuple(group_key.replaceFirst(/^.*.hrc\.r1-1\.ega\.grch37\./,""),file_list)}.set {ref_panel}

qc1_chr.filter{chrom, bfiles -> bfiles.collect {it.getBaseName()}.every {it.endsWith("step10")}}.map{group_key,file_list -> tuple (group_key.replaceFirst(/^*\.step10/,""), file_list)}.join(ref_panel).joi
n(genetic_map).into {shape_in;imp_in}
imp_in.view()
process phasing {
publishDir params.out, mode:'copy'
input:
tuple val(chrom), path(input_bed), path(input_ref), path(input_map) from shape_in
output:
tuple val("${chrom}"),path("${chrom}.*") into phased
script:
def (bed, bim, fam) = input_bed
def (haplotype, legend, sample) = input_ref
"""
shapeit --thread "${task.cpus}" --input-bed "${bed}" "${bim}" "${fam}" --input-map "${input_map}" --input-ref "${haplotype}" "${legend}" "${sample}" --duohmm --output-max "${chrom}.step10" --seed 54321
"""
}

```

# Answer
The merge operator is deprecated and will be removed from future Nextflow releases. You'll need instead the join operator to join by chromosome. You can use the fromFilePairs operator to quickly get the group key, formatting it as required:

```
channel
    .fromFilePairs("/some/path/HRC.r1-1.EGA.GRCh37.chr*.{hap.gz,legend.gz,samples}", size: 3)
    .map { group_key, file_list ->
        tuple( group_key.replaceFirst(/^HRC\.r1-1\.EGA\.GRCh37\./, ""), file_list )
    }
    .set { ref_panel_files }

channel
    .fromFilePairs("/some/path/genetic_map_chr*_combined_b37.txt", size: 1)
    .map { group_key, file_list ->
        tuple( group_key.replaceFirst(/^genetic_map_/, ""), file_list.first() )
    }
    .set { genetic_map_files }

qc1_chr
    .map { chrom, bfiles ->
        tuple( chrom, bfiles.findAll { it.baseName.endsWith("step10") } )
    }
    .join( ref_panel_files )
    .join( genetic_map_files )
    .set { phasing_inputs }

process phasing {

    input:
    tuple val(chrom), path(input_bed), path(input_ref), path(input_map) from phasing_inputs

    script:
    def (bed, bim, fam) = input_bed
    def (haplotypes, legend, sample) = input_ref

    """
    shapeit2 \\
        --thread "${task.cpus}" \\
        --input-bed "${bed}" "${bim}" "${fam}" \\
        --input-ref "${haplotypes}" "${legend}" "${sample}" \\
        --input-map "${input_map}" \\
        --duohmm \\
        --output-max "${chrom}.phased" \\
        --seed 54321
    """
}

```
> 3 votes

---
Tags: nextflow
---


# Question
Title: Removing duplicate FASTA sequences based on headers with Bash
I used the following command to remove duplicate FASTA sequences based on the header sequence:

```
paste -d $'\t' - - <fastaFileWithNoLinebreaksInSeq | sort -t $'\t' -uk1,1 | awk 'BEGIN{FS="\t";OFS="\n"}{print $1,$2}'

```

After this command, the output result is shown below:

```
>_10|_10
AAAGCAAAACACAGAATTGACGGAAAGACTTACGTTATTAAACGTGTTAAATATAATAACGAGTTTGGCATGGATTTTAAAGAAATAGAATTAATTGGCTCAGGTGGATTTGGCCAAGTTTTC
>_10|_11
ATTATGATCCTGAGACCAGTGATGATTCTCTTGAGAGCAGTGATTATGATCCTGAGAACAGCAAAAATAGTTCAAGGTTTGGCATGGATTTTAAAGAAATAGAATTAATTGGCTCAGGTGGATTTGGCCAAGTTTTCAAAGCAAAACACA
>_10|_12
TTTGGAACTCTTTGAACAAATAACAAAAGGGGTGGATTATATACATTCAAAAAAATTAATTCATAGAGATCTTAAGGTTTGGCATGGATTTTAAAGAAATAGAATTAATTGGCTCAGGTGGATTTGGCCAAGTTTTCAAAGCAAAACACA
>_10|_13
ACTTGTAACATCTCTGAAAAATGATGGAAAGCGAACAAGGAGTAAGGGAACTTTGCGATACATGAGCCCAGAACAGGTTTGGCATGGATTTTAAAGAAATAGAATTAATTGGCTCAGGTGGATTTGGCCAAGTTTTCAAAGCAAAACACA
>_10|_14
AACGAATTTCTTCGCAAGACTATGGAAAGGAAGTGGACCTCTACGCTTTGGGGCTAATTCTTGCTGAACTTCTTCATGTTTGGCATGGATTTTAAAGAAATAGAATTAATTGGCTCAGGTGGATTTGGCCAAGTTTTCAAAGCAAAACAC

```

This is my actual FASTA file: I don't want to lose header information from my FASTA sequence after removing duplicate sequences based on header information.

```
>ENST00000642616_2|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTAAAATGAAATCACAGGCGCATTACTGCCTTGTCTTGATGAGTCTCGTTTTGAAAATCTTGGAGTAAGCATGAAC
>ENST00000642616_2|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTAAAATGAAATCACAGGCGCATTACTGCCTTGTCTTGATGAGTCTCGTTTTGAAAATCTTGGAGTAAGGAAGTTG
>ENST00000642616_2|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTAAAATGAAATCACAGGCGCATTACTGCCTTGTCTTGATGAGTCTCGTTTTGAAAATCTTGGAGTAAGGGTGGGG
>ENST00000642616_2|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTAAAATGAAATCACAGGCGCATTACTGCCTTGTCTTGATGAGTCTCGTTTTGAAAATCTTGGAGTAAGGTAATTA
>ENST00000642616_2|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTAAAATGAAATCACAGGCGCATTACTGCCTTGTCTTGATGAGTCTCGTTTTGAAAATCTTGGAGTAAGGTCATGG
>ENST00000642616_2|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTAAAATGAAATCACAGGCGCATTACTGCCTTGTCTTGATGAGTCTCGTTTTGAAAATCTTGGAGTAAGTTCCTTG
>ENST00000642616_8|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTATTTGGTGCCTATCCTAAAACTTCCAGTGGGTTTTCTTGCCTTCCAGTGGCCATATAAAGGGCGTCCTGAAAAC
>ENST00000642616_7|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTCATGAACAAGGCTCAGTTATGATGTTTGAGCACCTTATTAATTCTAATGGAATTAAGCCTGTCATGGAACAATA
>ENST00000642616_1|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTGCGCCTGCGCGCGGGTCCGGCGCCGAGGTTCTTGACTGCTGTGCCGGACGCCAGGTGTAGCCATGCAGCGAGCC
>ENST00000642616_9|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTGGACTGCCATCATCTTGGAATCCAAAATAATTTTGATTACAAGCGCTTTATTAAGTTTGCCCGTGTCTGTGAAG
>ENST00000642616_5|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTGGTGGGGTATCTAGCAGGATGTCTAGTTCACGCACTGGGTGAAAAACAACCAGAGCTGCAGATAAGTGAACGAG
>ENST00000642616_4|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTGTAATTAATGATCCTATCCATGGCCACATTGAGCTCCACCCTCTCCTCGTCCGAATCATTGATACACCTCAATT
>ENST00000642616_6|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTGTCATGGGCCATTTTCTCACATGTTTGATGGACGATTTATTCCACTTGCTCGCCCGGAGGTGAAATGGACGATT
>ENST00000642616_6|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTGTCATGGGCCATTTTCTCACATGTTTGATGGACGATTTATTCCACTTGCTCGCCCGGAGGTGAAATGGACGCAT
>ENST00000642616_6|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTGTCATGGGCCATTTTCTCACATGTTTGATGGACGATTTATTCCACTTGCTCGCCCGGAGGTGAAATGGACGGAA
>ENST00000642616_3|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTTTCCTTGGGGGAGAGGAAGAAGCTGCTTAGTTATATCCAGCGATTGGTTCAAATCCACGTTGATACAATGAAGA
>ENST00000642616_3|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTTTCCTTGGGGGAGAGGAAGAAGCTGCTTAGTTATATCCAGCGATTGGTTCAAATCCACGTTGATACAATGAAGC
>ENST00000642616_3|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTTTCCTTGGGGGAGAGGAAGAAGCTGCTTAGTTATATCCAGCGATTGGTTCAAATCCACGTTGATACAATGAAGG
>ENST00000462861_2|ENST00000462861_3
AAAAAACTCTTCTACAGAAATTACTCTCAAAGAAACCTGAGGATCGACCTAATTTTTCACAGACCTACGGGATGGCATCATCTCAGATATATTTGATAAAAAAG

```

# Answer
If you want to avoid using extra libraries for any reason, you can just use a simple Python script (version 3.6 and above) to do this:

```
fr = open("dup_test.fasta", "r")
fw = open("dup_edited.fasta", "w")
seq_dict = {}
curr_header = ''

for line in fr:
   line = line.strip()
   if line[0] == '>':
        if line not in seq_dict:
          seq_dict[line] = ''
          curr_header = line
   else:
        if seq_dict[curr_header] == '':
          seq_dict[curr_header] = line
          fw.write(f"{curr_header}\n{line}\n")

fr.close()
fw.close()

```
> 10 votes

# Answer
You can use seqkit for this purpose.

`seqkit rmdup -n seqs.fa -o seqs_without_duplicate.fa`
> 9 votes

# Answer
One way using :

```
awk '/^>/ { f = !a[$0]++ } f' seqs.fa

```

---

**Explanation:**

The above is basically the same as:

```
awk '/^>/ { f = !($0 in a); a[$0]++ } f' seqs.fa

```

That is, for each header line (lines that start with a '\>' character), set a flag, 'f' to true if the current line is not in an array, 'a'. Then add the (header) line to the array. Note that in this form, it's not necessary to increment the value for each element - the following would be sufficient:

```
awk '/^>/ { f = !($0 in a); a[$0] } f' seqs.fa

```

The default action is to print, and we print lines only when 'f' is true. This is the same as:

```
awk '/^>/ { f = !($0 in a); a[$0] } f==1 { print $0 }' seqs.fa

```
> 4 votes

---
Tags: fasta, bash, linux
---


# Question
Title: Remove Redundant Sequences from FASTA file with Biopython reducing memory footprint
I need help with my Biopython code, I am kind of new to Python and Biopython so please don't be harsh on me.

Everything started from here: Remove duplicated sequences in FASTA with Python

I answered copying code from here: Remove Redundant Sequences from FASTA file in Python but I am still bugged about the OP's first question: Remove duplicated sequences in FASTA with Python

So I decided to try to minimize the code footprint on RAM. I ended up copying from Biopython documentation Split large file and producing:

```
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Nov  11 17:40:03 2020

@author: Pietro

"""

from memory_profiler import memory_usage

import memory_profiler

def batch_iterator(iterator, batch_size):
    """ from https://biopython.org/wiki/Split_large_file
    """
    entry = True  # Make sure we loop once
    while entry:
        batch = []
        while len(batch) < batch_size:
            try:
                entry = iterator.__next__()
            except StopIteration:
                entry = None
            if entry is None:
                # End of file
                break
            batch.append(entry)
        if batch:
            yield batch

# @profile
def f():

    from Bio import SeqIO

    import time

    import sys

    start = time.time()

    data_input_iter = SeqIO.parse("Newtextfile-03", "fasta")
    print('data_input_iter : ',sys.getsizeof(data_input_iter),'\n')

    c=0
    for i , batch in enumerate(batch_iterator(data_input_iter , 1)):
        record_iter = batch

        val = True
        c += 1
        if c == 1:
            with open("removed-double_Newtextfile-03", 'w+') as output_file:
                    SeqIO.write(record_iter, output_file, "fasta")
                    # print('copy   : ',(record_iter))
                    output_file.flush()
            output_file.close()

        else:

            SeqIO.write(record_iter, "temp_file", "fasta")

            record_iter = SeqIO.read("temp_file", "fasta")

            data_output_iter = SeqIO.parse("removed-double_Newtextfile-03", "fasta")

            d = 0

            try:
                for i , batch in enumerate(batch_iterator(data_output_iter  , 1)):
                        SeqIO.write(batch, "temp_file2", "fasta")

                        record_output = SeqIO.read("temp_file2", "fasta")

                        d += 1
                        if record_iter.seq != record_output.seq  :
                            # print('pass   : ',record_iter.name)
                            pass
                        else:
                            val = False
                            # print('false')
                            break

            except:
                print('ERROR !!!!!!!!!!!')
                quit()

            if val:
                        # print('copy   : ',record_iter.name)
                        with open("removed-double_Newtextfile-03", 'a+') as output_file:
                            output_file.write(str(record_iter.format("fasta")))
                        output_file.close()
            del data_output_iter

    end = time.time()

    print(f"Run time is {(end- start)}")

mem_usage = memory_usage(f)
# print('Memory usage (in chunks of .1 seconds): %s' % mem_usage)
print('Maximum memory usage: %s' % max(mem_usage))

# f()

```

but I have to major problems :

FIRST: I wasn't able to find a way to pass a SeqRecord object from a list to a file without saving it to a temp file first. That is without:

```
for i , batch in enumerate(batch_iterator(data_input_iter , 1)):
        record_iter = batch
......

SeqIO.write(record_iter, "temp_file", "fasta")

record_iter = SeqIO.read("temp_file", "fasta")

```

the line:

```
 if record_iter.seq != record_output.seq  :

```

will throw:

```
AttributeError: 'list' object has no attribute 'format'

AttributeError: 'list' object has no attribute 'seq'

```

you get first error and the if remove format from:

```
output_file.write(str(record_iter.format("fasta")))

```

and #quit() you get the second error How should I cope with this ? I mean is there a way to bypass the save and read lines ?

SECOND : and more important. I am sure my code doesn't use a set and a list .... and that's why is extremely slow compared to the copied Remove Redundant Sequences from FASTA file in Python posted below too. But how does it use even more RAM that the fastest code ? I tried to minimize the memory footprint but wasn't able !! Am I using a bad way to try to calculate the amount of RAM used ? And does Biopythom SeqIO.parse object already is the optimal choice for , lets say "2/4 Gb size" fasta files on a laptop with just 4Gb of RAM (remember the original OP post Remove duplicated sequences in FASTA with Python). Maybe my test file is not big enough ? Does anybody have a 2GB fasta file with a dozen of duplicates in it ?

Remove duplicated sequences in FASTA with Python memory usage:

```
size of parser :  48
memory seen :  32984
memory records :  2888
Run time is 0.0008047382036844889
Maximum memory usage: 55.359375

```

my code memory usage and speed:

```
data_input_iter :  48

Run time is 37.781240463256836
Maximum memory usage: 56.40234375

```

My fasta test set is: https://easyupload.io/z2yss1

Biopython method to remove duplicated from Remove duplicated sequences in FASTA with Python:

```
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
copied from https://bioinformatics.stackexchange.com/questions/13708/remove-redundant-sequences-from-fasta-file-in-python

uses set() instead of list for seen

should work faster for big files

"""

from memory_profiler import memory_usage

def f():

    import sys

    from Bio import SeqIO
    import time

    start = time.time()

    seen = set()
    records = []

    parser = SeqIO.parse("Newtextfile-03", "fasta")
    for record in parser:
        # print('running')
        if str(record.seq) not in seen:
            seen.add(str(record.seq))
            records.append(record)
            # print('copying')

    print('size of parser : ',sys.getsizeof(parser) )
    print('memory seen : ',sys.getsizeof(seen))
    print('memory records : ',sys.getsizeof(records))

    #writing to a fasta file
    SeqIO.write(records, "removed-double_Newtextfile-03-from-002", "fasta")
    end = time.time()

    print(f"Run time is {(end- start)/60}")

mem_usage = memory_usage(f)
# print('Memory usage (in chunks of .1 seconds): %s' % mem_usage)
print('Maximum memory usage: %s' % max(mem_usage))

# f()

```

# Answer
> 2 votes
I wrote the answer here and while I'm not clear on all aspects of this question, I can offer two suggestions to reduce the memory footprint considerably while maintaining high speed.

First, rather than appending the records to a list you could use a generator that generates elements on-the-fly rather than storing them in memory, e.g.

```
seen = set()

def yield_seqs():
    for record in SeqIO.parse("b4r2.fasta", "fasta"):
        if record.seq not in seen:
            seen.add(record.seq)
            yield record

SeqIO.write(yield_seqs(), "output.fasta", "fasta")

```

Second you could save even more memory, **but with the risk of false positives**, by using a bloom filter, e.g. using python-bloom-filter, just replace the set with `seen = BloomFilter(max_elements=10000, error_rate=0.001)`

# Answer
> 0 votes
FIRST: I wasn't able to find a way to pass a SeqRecord object from a list to a file without saving it to a temp file first. That is without:

```
for i , batch in enumerate(batch_iterator(data_input_iter , 1)):
        record_iter = batch
......

SeqIO.write(record_iter, "temp_file", "fasta")

record_iter = SeqIO.read("temp_file", "fasta")

```

the line:

```
 if record_iter.seq != record_output.seq  :

```

will throw:

```

AttributeError: 'list' object has no attribute 'seq'

```

After spending some timde on it I realized that:

```
print('\n@@@@@@@@@@@@@@\n',type(record_iter),'\n',record_iter)

```

is returning:

```
<class 'list'>

```

but jus using :

```
for i , batch in enumerate(batch_iterator(data_input_iter , 1)):
        record_iter = batch[0]

```

now :

```
print('\n@@@@@@@@@@@@@@\n',type(record_iter),'\n',record_iter)

```

is returning:

```
<class 'Bio.SeqRecord.SeqRecord'>

```

no need to save and read the Bio.SeqRecord.SeqRecord' anymore

so this is just half the question answered so no need to save and read the record on file

but even if speed increases the memory usage goes up :

data_input_iter :  48

Run time is 1.442115306854248
Maximum memory usage: 110.09765625
```

---
Tags: python, fasta, seqio, biopython
---


title: How to convert txt.gz into fastq.gz?

text:
# Question
Title: How to convert txt.gz into fastq.gz?
I'm a newbie to bioinfo. I'm trying to use galaxy.networkanalyst.ca for the analysis of RNA-Seq data and for that I need to convert ***txt.gz files into fastq.gz files***.

Link: GSM4273445 is one sample and the whole series is from here.

I've Googled to find out if there are online tools to convert it, but there aren't any. Besides, there is no discussion addressing this format conversion on researchgate.

I intend to use this RNA-Seq data and investigate the relative expression of some specific transcripts of interest. I have been introduced to 'galaxy.networkanalyst.ca' where the very first step is to upload fastq files. In addition to this I was aware of ncbi GEOdatasets as source of open RnA-Seq data.

Can anyone suggest a method or tool for this?

# Answer
> 3 votes
I've had a quick look at the file you linked to and this is a featureCounts file - this list the number of reads for each gene/transcript (COUTNS.txt in the file name pretty much gave it away that it was going to be a file of, well counts).

I can't help you with Galaxy, but if you want to do an RNA-seq analysis, the count data + metadata is likely all you'll need.

If you really want to do the preprocessing and mapping yourself, the fastq files can be downloaded from the SRA archive (SRP242087) - click on the SRA link in your whole series link.

---
Tags: fastq
---


title: Error downloading GSE dataset

text:
# Question
Title: Error downloading GSE dataset
I've installed GEOquery with

```
install.packages("BiocManager")
BiocManager::install("GEOquery")
library(GEOquery)

```

then I ran the command

```
GSE68808 <- getGEO("GSE68808")

```

but it resulted in..

```
Found 1 file(s)
GSE68808_series_matrix.txt.gz
trying URL 'https://ftp.ncbi.nlm.nih.gov/geo/series/GSE68nnn/GSE68808/matrix/GSE68808_series_matrix.txt.gz'
Content type 'application/x-gzip' length 4138819 bytes (3.9 MB)
downloaded 3.9 MB

error: 'C:\Users\(myname)\AppData\Local\Temp\RtmpUp5Sv7\file8e844a7a59bc' does not exist.

```

ㅜㅜ

How can I solve this? thanks.

# Answer
> 0 votes
There seems to be some problem with your "platform specific tempdir" that is used by getGEO's `destdir` parameter. See `?getGEO` for more info on the `destdir` parameter.

Try

```
GSE68808 <- getGEO("GSE68808", destdir = getwd())

```

---
Tags: geoquery, gse
---


# Question
Title: Is there a good practice (or an easy one for a SBOL newbie) to set up a collection of biological parts and devices?
I have a collection of biological parts and devices in a particular format that basically stores information of the sequence of the elements. My intention is to move to a more standardized data structure with a particular interest in Synthetic Biology. SBOL 3 looks like a good choice, then the question is if there is a way to kind of "import" data from my sequence-based database into SBOL?

# Answer
For SBOL 2, you can readily import data via the online converter tool at https://converter.sbolstandard.org/ The converter basically just wraps the conversion function found in the SBOL libraries, so if you're coding, you can do the conversion directly from the libraries.

For SBOL 3, there isn't yet an online tool, but James McLaughlin has a converter in his JavaScript library: https://github.com/sboltools/sbolgraph I believe he is working to bring it online shortly, so hopefully this answer will soon need to be updated or replaced by an answer from him.

An anecdote on the difficulty level: I recently started using an SBOL3--\>SBOL2--\>GenBank conversion using pySBOL3 and pySBOL3. It took a couple of hours to setup the pipeline, and isn't perfect, but it's pretty doable (and ought to be generalized for others to use as well at some point). The other way would be easier, since GenBank contents are pretty simple.
> 5 votes

# Answer
Further to jakebeal's answer there is also an excel2sbol library (still in alpha) that can take input to a standard template and convert to SBOL2.

For SBOL3 there is also a python library, but using this to create SBOL3 is more involved.
> 3 votes

---
Tags: sbol
---


# Question
Title: How to use SBOL (Synthetic Biology Open Language)?
One of my students is working on the development of a SynBio tool to design and represent genetic circuits in a standardized format. I know there are some libraries out there.

What I am interested in finding out is to have a hands-on tutorial how to use an SBOL library (say, pySBOL) to represent genetic circuits in SBOL data and visual representations.

# Answer
For pySBOL there is a tutorial you can follow in the documentation. Additionally, there is a youtube video walking through a similar tutorial, note the tutorial starts at about 1:09:00.

There are lots of additional tools and libraries out there such as SBOL Canvas (for visualisation), libsbolj (a java library for SBOL), and other libraries and tools listed on the SBOL website.
> 3 votes

---
Tags: sbol
---


# Question
Title: How to trigger next process on completion of previous process?
Simple question but I am stumped. (I am still learning Nextflow)

I have a process which takes as input a channel of files and processes each file separately.

There is a second process which needs to wait until all files are complete, so it can merge the outputs into a single collated file.

How can I do this? I have tried subscribing to "onComplete" event to no avail as I cannot specify a process as the function to trigger.

files = Channel.from params.the_files_list

process process_files {

    input:
    val file_path from files

    output:
    stdout result

    """
    python ./process_one_file.py $file_path
    """
}

process merge_results {

    input:
    ????

    output:
    ????

    script:
    """
    echo "How to trigger this process once, after all individual files are processed?"
    """

}

??? Want to do something like this
result.subscribe onComplete: merge_results

```
</code></pre>

# Answer
> 2 votes
So from this page... https://groups.google.com/g/nextflow/c/jt77\_-uApMs

It is the "scatter-gather" pattern. I think the answer is to collect() the channel as follows - interested to hear if there is a better way

```
files = Channel.from params.the_files_list

process process_files {

    input:
    val file_path from files

    output:
    stdout result

    """
    python ./process_one_file.py $file_path
    """
}

process merge_results {

    input:
    val go from result.collect()

    script:
    """
    echo "Triggered once after all files complete!"
    """

}

```

---
Tags: nextflow, workflow-management
---


# Question
Title: About optimal distance (or dissimilarity) measure in proportion data
My data is proportional data, but variables are correlated with each other. So I am wondering what kinds of distance or dissimilarity we should use to take correlation into account.

If we would like to use data with proportion (e.g.,total sum scaling or compositional data), what kinds of distance or dissimilarity can we use? Bray-Curtis is one of them, but does this make sense to use arcisin or logit transformation and use Pearson correlation?

In continuous variables, we can use Pearson and Spearman distance. In our case, what kinds of distance can we use? Pearson distance cannot be applied to data with 0-1 scale. Logit transformation -\> Pearson?

# Answer
I would suggest reading a little on the topic; my understanding is that distance measures can be pretty robust to proportion/compositional data. Here is a relevant textbook chapter:

> Other distance measures assume that the data are proportions ranging from zero to one, inclusive.

Among these are indeed Bray-Curtis as well as Jaccard distance (and others). I think that you will want to use one of the more traditional distances rather than a Pearson/Spearman, though they are getting at the same problem. A transformation might help, I suppose, but I would prefer to choose a simple but predictable distance measure rather than a correction on a possibly unsuitable distance that might have unpredictable behavior.

It's a little hard to tell from your question, but it seems that you are wanting to in some way account for indirect correlations using a distance metric. I would not suggest trying to make your distance measure do all this work. Instead, I would use a simple but clear distance method, and then use a multidimensional analysis such as principal coordinate analysis (or similar) to understand the correlation structure of your compositional variables. This is a fairly standard workflow in metagenomics, for instance, which relies heavily on high-dimensional compositional data. I would recommend looking more into metagenomics statistical applications for inspiration if needed.

The choice does also depend to a certain extent on the specific application. In some fields certain methods are familiar and trusted even if there is not a strong empirical/theoretical rule.

You've asked a lot of questions about Pearson/Spearman correlations. Much has been written on the use of correlation analysis on compositional data, most of it cautioning against doing this kind of analysis without special treatment. This leads to difficulties in building meaningful association models for different variables. One method that I'm familiar with in metagenomics is sparse correlation, e.g. SPARCC. But there are newer methods that may be better, for example here and here.
> 1 votes

---
Tags: r
---


# Question
Title: Extract residue sequence from pdb file (in biopython, but open to recommendations)
I'm new to Biopython and I'd like to extract the sequence of residues from a pdb file.

My two questions are:

1. What is the simplest way to do this? (Esp. when there is more than one sequence) and
2. Should I be worried about "PDBConstructionWarning: WARNING: Chain B is discontinuous"?

So far, I've obtained the residue sequence via:

```
p = PDBParser()
structure = p.get_structure("1ppi", "1ppi.pdb")
ppb=PPBuilder()
for pp in ppb.build_peptides(structure):
     print(pp.get_sequence())

seq = pp.get_sequence().__str__()

```

This seems to work well for this molecule. However, is there an easier way, especially when there is more than one sequence?

For example, I've read that one can also do

```
res_list = Bio.PDB.Selection.unfold_entities(structure, 'R')

```

but res\_list is not a sequence of residues in str, and I don't know how to convert the output from res\_list to a str sequence.

In addition (or perhaps because I'm going through the PPBuilder), I've recently gotten a lot of warnings of the type: /usr/local/lib/python3.8/site-packages/Bio/PDB/StructureBuilder.py:89: PDBConstructionWarning: WARNING: Chain A is discontinuous at line.. (For example, with `structure = p.get_structure("5owu", "5owu.pdb")`)

I've seen a discussion about how to silence such warnings, but should I be worried about this? I've also noticed that when I get these warnings, pp builder seems to give me more sequences then are there.

# Answer
> 4 votes
1. There are other ways to do it (many and some might be easier), but if you want to do it with BioPython PDB module you can iterate the residues like this:

```
 # You can use a dict to convert three letter code to one letter code
 d3to1 = {'CYS': 'C', 'ASP': 'D', 'SER': 'S', 'GLN': 'Q', 'LYS': 'K',
 'ILE': 'I', 'PRO': 'P', 'THR': 'T', 'PHE': 'F', 'ASN': 'N',
 'GLY': 'G', 'HIS': 'H', 'LEU': 'L', 'ARG': 'R', 'TRP': 'W',
 'ALA': 'A', 'VAL':'V', 'GLU': 'E', 'TYR': 'Y', 'MET': 'M'}

 # Just an example input pdb
 record = '1pa2.pdb'

 # run parser
 parser = PDBParser(QUIET=True)
 structure = parser.get_structure('struct', record)

 # iterate each model, chain, and residue
 # printing out the sequence for each chain

 for model in structure:
     for chain in model:
         seq = []
         for residue in chain:
             seq.append(d3to1[residue.resname])
         print('>some_header\n',''.join(seq))

```
2. Not really, is very common, you can turn it off like this (taken from https://lists.open-bio.org/pipermail/biopython/2014-July/015371.html):

```
 import warnings
 from Bio import PDBConstructionWarning
 #your code
 with warnings.catch_warnings():
 warnings.simplefilter('ignore', PDBConstructionWarning)
 #your code which might trigger the warning
 #rest of your code here

 However, because this is so common, you can just use the
 QUIET=True option on the PDBParser object:

 from Bio.PDB.PDBParser import PDBParser
 struct = PDBParser(QUIET=True).get_structure('tmp', pdbf)

```

Try help(PDBParser) for more, or see: http://biopython.org/DIST/docs/api/Bio.PDB.PDBParser%27.PDBParser-class.html#\__init\__

# Answer
> 7 votes
The discontinuity is normal: it is a missing stretch due to lack of density. In PyMOL you see them as dotted lines.

I am sorry to say, the easiest/safest way to get the sequences is not using the PDB files...

```
import requests
data = requests.get(f'https://www.ebi.ac.uk/pdbe/api/pdb/entry/molecules/{code}').json()[code.lower()]
print(data[0]['sequence'])

```

One gets back a lot of information: example. This will tell you when you have multiple copies of the same peptide. The only information you do not get is what the PDB numbering is relative to the canonical Uniprot for that chain, which has to be taken from SIFTS.

A side note... Biopython.PDB is one of the best BioPython packages, but I personally I much prefer PyMOL as a python module (`pymol2` module) to Biopython.PDB as it is functionally more complete and straightforward even if everything operates in a stage. You might want to check it out.

```
with pymol2.PyMOL() as pymol:
    pymol.cmd.fetch('1UBQ', 'prot')
    print(pymol.cmd.get_fastastr('prot'))

```

# Answer
> 4 votes
You can also use the SeqIO module from Biopython https://biopython.org/docs/1.75/api/Bio.SeqIO.PdbIO.html

and do something like this:

```
from Bio import SeqIO

for record in SeqIO.parse(target, "pdb-atom"):
    print(record.seq)

```

where target is the path to your pdb file the parser has multiple input/output functions, and pdb-atom is the one that will read out of the coordinates itself.

# Answer
> 0 votes
you can try PDBx Python Parser too:

http://mmcif.wwpdb.org/docs/sw-examples/python/html/index.html

it has an example showing exactly what you need

http://mmcif.wwpdb.org/docs/sw-examples/python/html/fasta.html

---
Tags: python, pdb
---


# Question
Title: CRISPR/Cas9 screen analysis with Mageck: paired-end sequencing
I want to analyse data from a CRISPR/Cas9 screen (control vs. treatment) and I'm using Mageck (https://sourceforge.net/projects/mageck/). The problem is that I'm working with paired-end sequencing data, do I have to analyse R1 and R2 fastq files separately or is possible to do it at same time with Mageck?

# Answer
The other answer is incorrect. The original poster is asking about paired-end sequencing data, not paired treatment/control experimental design.

Here is how you analyze paired-end sequencing data (this was an answer I gave on Biostars last year):

https://www.biostars.org/p/303656/
> 1 votes

---
Tags: ngs, crispr
---


# Question
Title: Manta SV caller: I can't find a file in the path specified by the user guide
I trying to locally reproduce a Manta run as seen in the *log* of an illumina basespace RNA-Seq Alignment tool v1.1.1 run I've done. The log reads as follows:

```
[LOG][2021-03-12 20:58:56.377]03/12/2021,20:58:56.193,Launch process:  \
/path_to_RNA_Tools/bin/python  \
'/path_to_RNA_Tools/Manta/bin/configManta.py'  \
--rna --normalBam='path_to_alignments/alignments.bam' \
--runDir='path_to_samples/CTR/replicates/CTR/MantaFusions/manta.run'  \
--referenceFasta='/path_to_genomes/genome.fa' \
--config '/opt/illumina/Isis/2.6.25.18/Workflows/WholeGenomeRnaSeqWorker/scripts/configManta.py.ini'

```

The file shown in the line above `configManta.py.ini` is however no were to be found. The Manta user guide states in the linked section

> Options listed in the file: `${MANTA_INSTALL_PATH}/bin/configManta.py.ini`

The `configManta.py` is in the expected directory, only the `configManta.py.ini` is absent.

I am using Manta version 1.6.0 from bioconda.

Your insights will be much appreciated!

# Answer
> 2 votes
The `--config` file option for Manta can be used to override the defaults in the global config file. If you type `configManta.py -h`, you'll see where the global config file has been installed. For me, it's here:

```
/home/steve/miniconda3/envs/manta/share/manta-1.6.0-1/bin/configManta.py.ini

```

So the command in your log file takes a specific config file installed into:

```
/opt/illumina/Isis/2.6.25.18/...

```

I suspect 'Isis' is a reference to Illumina's analysis tool 'IsisRNA'. Papers that mention Illumina's 'RnaReadCounter' tool, say that it is part of the *internal* analysis tool 'IsisRNA'. Because it's one of Illumina's in-house tools, I think you'd need to contact Illumina support to get a copy. I would only do this though if you need to reproduce *exactly* as per the command in your log file. If you only need to roughly reproduce the output, the defaults in the global config file are probably fine.

---
Tags: fusions
---


# Question
Title: Inflate operation failed: progress temporarily not possible, or in() / out() returned an error
I download a BAM file from ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics\_ChromiumGenome\_LongRanger2.0\_06202016/NA12878\_GRCh38.bam (171G), and I am trying to view the comparison result of the BAM file using the following command:

```
[yyn@localhost All]$ samtools flagstat NA12878_GRCh38.bam

```

The result shows that the file is truncated, the result as follows:

```
[E::bgzf_uncompress] Inflate operation failed: progress temporarily not possible, or in() / out() returned an error
[E::bgzf_read] Read block operation failed with error 1 after 0 of 4 bytes
[bam_flagstat_core] Truncated file? Continue anyway.
72639 + 0 in total (QC-passed reads + QC-failed reads)
896 + 0 secondary
0 + 0 supplementary
2529 + 0 duplicates
72639 + 0 mapped (100.00% : N/A)
71743 + 0 paired in sequencing
36865 + 0 read1
34878 + 0 read2
48600 + 0 properly paired (67.74% : N/A)
69949 + 0 with itself and mate mapped
1794 + 0 singletons (2.50% : N/A)
10233 + 0 with mate mapped to a different chr
2869 + 0 with mate mapped to a different chr (mapQ>=5)

```

The total number of reads is only 72639, Obviously, this is the wrong result. And I checked the file again with the following command

```
[yyn@localhost All]$ samtools quickcheck -qv NA12878_GRCh38.bam && echo 'all ok' || echo 'failed check'

```

The result is:

```
all ok

```

I checked the contents of the file again, and the result showed that there was no problem with the first 72620 lines in the file, but there was a problem with the first 72630 lines, so the file should be truncated in the 72620 to 72630 lines.

```
[yyn@localhost All]$ samtools view NA12878_GRCh38.bam | head -n 72610 >> head72610.sam

```

That's ok.

```
[yyn@localhost All]$ samtools view NA12878_GRCh38.bam | head -n 72620 >> head72620.sam
[E::bgzf_uncompress] Inflate operation failed: progress temporarily not possible, or in() / out() returned an error
[E::bgzf_read] Read block operation failed with error 1 after 0 of 4 bytes
[main_samview] truncated file.

```

`samtools sort` and `samtools index` are the same error results as above. Does anyone know what's going on? Thanks!

# Answer
My best guess would be a corruption during download as:

```
samtools quickcheck -qvvvv ftp://ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/NA12878_GRCh38.bam
verbosity set to 4
checking ftp://ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/NA12878_GRCh38.bam
opened ftp://ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/NA12878_GRCh38.bam
ftp://ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/NA12878_GRCh38.bam is sequence data
ftp://ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/NA12878_GRCh38.bam has 196 targets in header.
ftp://ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/NA12878_GRCh38.bam has good EOF block.

```

Plus I can successfully run `idxstats` from remote and fetch individual chromosomes without error via `samtools view`. Try redownloading before digging deeper.
> 1 votes

---
Tags: sam, samtools
---


# Question
Title: Why is there no GT-entry in my .vcf file from bacteria
I mapped raw reads obtained from an E. coli genome (.fastq file) to my reference genome. Next, I exported the .vcf (variant calling file), which looks something like the header example shown below (total \>50 lines):

```
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Variants:_LGC19-XL01_S63_L001_R_001_(trimmed)
CFC381_K12_Bw25113  360287  .   GGGAAT  AATTCC  1952.76 .   NS=1;VF=0.380;SB=0.533;SB50=0.60;SB65=1.0;TYPE=Substitution;AVQUAL=30   DP:AO   231:89

```

**How can I back-fill or quickly assign the missing GTs to the `.vcf` ?**

For my downstream process (`plink --pca`, via `bfctools merge`), I need that there are GT entries in every `.vcf`, so I can merge them and the respective GT is shown in the merged `.vcf`. I only found solutions that need the `.bam` file, but that is laborious and I don't think I need the `.bam` in the case of haploid data.

However, the`.fastq` is actually from a mixed culture, so I wonder do I really need to follow the `.bam`-strategy? Suggestion for a tool? But in that case, why did it not export the GTs in the fist place (I exported from "geneious" sequence analysis program)

# Answer
> 2 votes
The representation for a variant present in a haploid genome should be `1`. A heterozygote in a diploid genome would be `0/1` (or `0|1` or `1|0` if phased). So, to add a `1` as the genotype, you would need to first insert a `GT` to the beginning of each INFO line and then the `1` at the beginning of the last field. But, this should happen only on lines that are not headers (do not start with `#`). This should do it:

```
$ awk -F'\t' -v OFS="\t" '{ if(/^#/){print}else{$9="GT:"$9;$10="1:"$10;print}}' file.vcf
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Variants:_LGC19-XL01_S63_L001_R_001_(trimmed)
CFC381_K12_Bw25113  360287  .   GGGAAT  AATTCC  1952.76 .   NS=1;VF=0.380;SB=0.533;SB50=0.60;SB65=1.0;TYPE=Substitution;AVQUAL=30   GT:DP:AO    1:231:89

```

So just redirect the output to a fixed file:

```
$ awk -F'\t' -v OFS="\t" '{ if(/^#/){print}else{$9="GT:"$9;$10="1:"$10;print}}' file.vcf > fixed.vcf

```

Now, in this specific case, you could also do something like this:

```
$ sed -i.bak '/^[^#]/{s/\<DP:AO\t/GT:DP:AO\t1:/}' file.vcf

```

But that assumes that all lines will have the exact same INFO field of `DP:AO` which is not a safe assumption since each variant can have different INFO fields.

---
Tags: variant-calling
---


# Question
Title: scRNA-seq multi-dataset integration for small datasets
There have been a few methods proposed for integration (or batch correction) of scRNA-seq datasets, such as Seurat CCA, MNN Correct, Scanorama, and Harmony. The concern is generally about the maximum number of cells that they handle, but I haven't seen any discussion about the minimum number of cells. I am confident they can all handle 10k cells reasonably well and will fail with 10 cells, but where do you draw the line? Is there a method that works best for small datasets?

For example, with plate-based platforms like Fluidigm, many experiments only have 96 cells and potentially much less after quality filtering. How can those be used?

# Answer
> 1 votes
TLDR: Harmony can work on 106~ cell samples but has less versatility then methods like BATMAN. BATMAN is useful if you need your data for differential gene expression and single-cell eQTL and can work on ~200 cell data.

Well their certainly is a minimum which comes down to several factors. As many of these rely on algorithms like KNN and other clustering algorithms for integration the threshold is dependent on the quality of the data and the algorithms abilities to predict reasonable clusters from that data. As well as the downstream uses of the data...

BATMAN highest performance on simulated set 200 cells and 1000 cells
**With 200 cells from one batch and and 1000 in the other only BATMAN succeeded in clustering the high dimensionality simulated data well. If you need your data for differential gene expression and single-cell eQTL, BATMAN seems to be the correct choice, for small sets.** "The downside of these methods is that they operate in latent space, which limits their interpretability and use in downstream analyses such as differential gene expression and single-cell eQTL analyses." "A shows that the two original datasets become better mixed with each other. BATMAN has the highest performance not only when considering the top two principal components but also when considering more top principal components. It is the only method that manages to efficiently maintain a high iLISI score in a larger number of dimensions."

Harmony is best if... differential gene expression and single-cell eQTL is not necessary **"Moreover, we show that Harmony requires dramatically fewer computational resources. It is the only available algorithm that makes the integration of ~106 cells feasible on a personal computer."**

Computing time benchmarks-section It seems as the runtime of these algorithms is less than 2 minutes, using a set with less than ~2000 cells will not benefit your computational time much. "When processing a small data set of ~ 2000 cells, all four methods took less than 2 min." "To obtain such datasets, we downsampled the MCA and TM datasets to obtain a total of 9 sets of data containing between ~ 2000 and ~ 140,000 cells, while the number of highly variable genes (HVGs) was controlled in a range from ~ 2000 to ~ 3000 (Table S1)."

---
Tags: scrnaseq
---