# Question
Title: What is the most appropriate way to find the most recent common ancestor between two distantly related species
I want to specifically find the common ancestor between a lobster and a humans. I suspect it was an aquatic worm of some description. But I want to know about the nervous system of this common ancestor.

# Answer
> 4 votes
Because I've now posted several comments, I'll just roll them all up.

For background on the approaches used to identify most recent common ancestors and a high-level look at how animal taxonomy has been inferred, I suggest Lynch 1999.

I think that there are 2 interpretations of this question. If you are interested in **just looking up a single MRCA of well-defined clades**, such as lobster and human, here are some approaches:

## Easy way:

1. Look at a tree diagram, e.g. this:

2. Find the tips that correspond to your species of interest (arthropods for lobster, chordata for humans).
3. Find where they join together in the diagram (the branch labeled "true coelom").

You have your answer, the MRCA is the group of organisms with a true coelom, coelomates.

## A more involved way using a database

1. Go to this website.
2. Find the group of species 1 (arthropods, protostomes, etc. for lobster, chordata, deuterostomes etc. for human)
3. navigate around until you see the group containing the two groups (in this case listed as "bilateria"). In this case you are looking for the bilaterian common ancestor.

## another database

1. Go to this website.
2. Point and click your way to a view where you see your 2 clades of interest (arthropods, chordates in this case). See figure.
3. Find where they join (in this case, it is less certain about the existence of a coelomate common ancestor, so it just says "**bilaterians**").

## In R

We might also want a programmatic procedure in case we need to automate it. The following gives a series of commands in R that do this.

Get a newick/nexus tree from somewhere e.g. here of animals (go "control panel --\> export, format = "newick tree" --\> download).

For this tree, we can approximate an answer as following:

```
library(ape)
it = read.tree("~/Downloads/zUtVMYn3wJQBoP7E8fygTw_newick.txt")  # autogenerated filename from the website.
getMRCA(it, c("Homo_sapiens", "Drosophila_melanogaster"))  # drosophila is one of only 2 arthropods in this tree, but arthropods all have an arthropod common ancestor
# prints "200", the node index of the MRCA.

```

You can use that node index to do things like look up a name in the `node.label` attribute of the tree (`it`) structure, or to do other operations on the tree.

# Answer
> 3 votes
Good answer.

The common ancestor is not known due to the paucity of the fossil record at this point in evolutionary time. A reconstruction can be made using coalescent theory which is what @MaximillanPress describes. The level of uncertainty at that depth of calculation requires significant expertise. I would re-evaluate the comparison with related organisms that are heavily studied.

Lobsters are decapods (Decapoda) and closely related to insects, moreover they are monophyletic with the Order Insecta. A lot is known about a lot of insects, many being model organisms, such as Drosophila melanogaster. Both decapods and insects exist within the arthropods (Arthropoda), along with spiders and the like. Their ancestors were marine arthropods, trilobites are arthropods. The major diversity of the arthropods was generated in the Cambrian Explosion and a lot of work has gone into identifying these fossils.

Traversing between decapods and insects is as far as I'd go, which is likely the late stages of the Cambrian explosion. There are likely to be putative ancestral fossils - not that you can get any molecular sequence from them however. The closest living ancestor on the insect side is Colembola (beach hoppers).

There are serious difficulties in major molecular calculations between decapods and mammals.

---
Tags: sequence-alignment, phylogenetics, genome, phylogeny
---


# Question
Title: What tool RCSB and AlphaFold use to visualize 3D structures?
I would like to implement on a web page a tool to visualize 3D structures like it is for AlphaFold and RCSB. A couple of examples: https://www.rcsb.org/3d-view/7KX9 https://alphafold.com/entry/Q5VSL9 As you can see, if you hover over a part of the structure or on a residue on the RNA/amino-acid sequence, that residue is underlined and different details pop up on the screen.

What is the javascript library that they use? Do you know also where to find the documentation? Best regards!

# Answer
Different criteria give different rankings.

* `Mol*` (read molstar with trilled rhotic R according to the given IPA) is the newest, is used by the RCSB PDB and can support huge complexes. It is less implemented and has a tricker documentation.
* `NGL` is the former viewer from RCSB PDB. It is good and the switch to Mol∗ was driven by an effort for uniformity with the PDBe. Both are written by Alex Rose. There is more documentation about NGL —in the form of a gallery of codepens and various GitHub repos using it. There is a package for Python notebooks called nglview. It has issues such as no SS predictions (i.e. requires sheets and helix entries to be defined in the file for proper rendering).
* `3DMol` is simpler and is of note as the Python module nglview does not work in colab notebooks, but the Python module py3DMol does.
> 5 votes

---
Tags: protein-structure, 3d-structure, rna-structure
---


# Question
Title: Analysis of differential transcript usage (DTU)
Recent breakthroughs in bioinformatics tools for quantification (e.g. Cufflinks/Kallisto/Salmon etc.) and tools which can identify differential transcript usage (DTU) (e.g. DRIMSeq, Cufflinks etc.) mean that from RNA-seq data we can now relatively easy obtain a genome wide analysis of transcripts that are differentially used betwen conditions.

What can you use these results for? What systematic analysis does this transcript level data enable?

# Answer
> 4 votes
I will take the liberty of giving one possible answers to my own question – but I’m very interested in other answers.

One analysis type that such data enables is the analysis of transcript switches with predicted potential consequences.

I myself have recently developed such a tool called IsoformSwitchAnalyzeR. IsoformSwitchAnalyzeR enables statistical identification (via DRIMSeq) of isoform switches with predicted functional consequences. The consequences analyzed can be chosen from a long list which includes gain/loss of protein domains, signal peptides changes in NMD sensitivity etc. The R package also enables easy visualization of isoform switches along with their consequences and it directly supports the output of Cufflinks/Cuffdiff, RSEM, Salmon and Kallisto.

Apart from enabling identification of interesting examples switch identification also enables systematic analysis of what genes are affected. For inspiration I will recommend one of my own articles (describing results obtained with IsoformSwitchAnalyzeR) as well Hector et al’s recent paper (which is not using IsoformSwitchAnalyzeR). Of particular interest and finesse is Hector’s analysis of how isoform switches can disrupt protein-protein interactions.

Looking forward to hear more ideas.

# Answer
> 3 votes
I am not sure if this has been done, is common, or are research lines, but here is what I think can be done with transcripts' differences (aside from comparing the change of the expression of each transcript):

Having the different transcripts, enables us to study also the relation between transcripts, which are co-expressed together or differently co-expressed. Maybe at the individual level the switch between transcripts of a gene is not related to a condition, but together with other transcripts it is.

Transcripts usage could be also linked to structural or genomic changes. Either short/small changes or big. At genomic level SNP, inversions or methylation..., or structural changes such as compression of a region, different distances between chromosomes/regions...

# Answer
> 1 votes
Additional insights can frequently be made in research by studying things at a lower level. By peeling back the higher layers, effects can be seen that are hidden (or reversed) when looking at aggregate data. This is generally referred to as Simpson's Paradox.

We know that RNA is functional at an isoform/transcript level, and in some cases the expression of different transcripts from the same gene can have different biologically-relevant outcomes. By looking at transcript-level information, even if the transcripts for a gene average out to be a similar level between different treatment conditions, variation at an individual transcript level can still be discovered. This is important when looking at differential expression of *populations* of cells, because it can clue researchers into finding out that a particular cell sub-population might be of interest; single cells typically will only produce single isoforms for each gene at a specific point in time.

# Answer
> 0 votes
I understand that this post is rather old but there has been some development in the field and this 2019 paper might be a good resource for interested individuals.

---
Tags: rna-seq, transcriptome, software-recommendation
---


title: 'GCLOUD_SDK_PATH' Error

text:
# Question
Title: 'GCLOUD_SDK_PATH' Error
I have run into the following problem

```
Error: failed to find 'gcloud' binary; set option or environment variable 'GCLOUD_SDK_PATH'

```

when executing this line of R code:

```
avfiles_ls( recursive = FALSE,
            namespace = avworkspace_namespace('bdcat-cohort3-stanislawski'),
            name = avworkspace_name('SOLGeneticTopMed')
)

```

This works when I run it in a notebook in a cloud environment but not when I am running an R script inside a wdl on Terra.

Is there a way to more fully quality where the namespace is?

Thanks. Liz

# Answer
> 2 votes
I assume that you are running Terra on Google Cloud Platform (GCP)?

It is weird, because Terra is supposed to run on GCP - that is what its intended purpose - and finding gcloud is automatic, you don't need to configure it.

Three suggestions:

1. If the cloud is GCP it is appears to could be a bug with the Broad Institutes platform and should be reported directly on their website. They will know immediately, so if they don't respond immediately its a user problem. This is because they will prioritise system bugs, but user errors they may not respond to.
2. At a guess the error suggests you just need to point to gcloud within .bashrc or .profile, but you don't need to do this on GCP because gcloud is its modus operandi\*.
3. Run WDL directly outside Terra.

*The simplest approach is to install miniconda and use it to install gcloud-sdk package and then install Terra from within conda. Assuming Terra is present on "-c bioconda". At this point Terra is forced to see gcloud*.

However, it feels like Terra is being installed on a cloud outside GCP and hence the error.

* To those who don't know GCP gcloud = google cloud.

---
Tags: python
---


title: How to make MSA for contact prediction?

text:
# Question
Title: How to make MSA for contact prediction?
I'd like to try a contact prediction method (like TripletRes, CCMpred etc) using a custom set of sequences. There are about 1.5M sequences, all very similar to the canonical sequence; their distribution reflects the occurence of mutations "in the wild" (they're not clustered or anything like that). How should I make a MSA in this case?

For example in the original **TripletRes paper** it looks like the MSA comes from DeepMSA - this searches UniClust30 and UniRef90 using **hhblits** and **jackhmmer**. I don't see any way to use DeepMSA to make a MSA from a fully custom set of sequences, and also it's not clear how to combine a custom set with the UniClust/UniRef search results. I can make an alignment of just my data with jackhmmer, but I don't know if that's what the other tools expect. For one thing, the alignment is pretty huge compared to the typical inputs (~million vs typical DeepMSA outputs of a few hundred); for another, it is a lot less variable (most sequences differ in just 1-2 aa, vs completely different species).

# Answer
> 1 votes
That's a LOT of sequences! I'm doubtful you'd be able to get TripletRes or CCMPred to run at all with an MSA depth of 1.5 million. I suggest pruning your dataset to remove redundant sequences (sequences which don't add much information, i.e. the ones which differ only 1 or 2 AA). I'm curious if your Neff (normalized number of effective sequences--as defined in the ResTriplet paper you link) is anywhere close to your number of sequences.

Also, I think you've misunderstood the use of DeepMSA. DeepMSA's main use is to find sequences to include in an MSA. I believe it uses much simpler programs to construct the MSA itself. So you'd use DeepMSA if you wanted to do 1 query sequence --\> an MSA of several hundred sequences. Since you already have the (very complete) set of sequences, it wouldn't make sense to use DeepMSA. You could instead just use something like Clustal Omega to generate the MSA. There's also a pretty easy-to-use web server but it's limited to 4k sequences.

# Answer
> 1 votes
You might have some luck with Parsnp. From the description:

> Parsnp was designed to align the core genome of hundreds to thousands of bacterial genomes within a few minutes to few hours. Input can be both draft assemblies and finished genomes, and output includes variant (SNP) calls, core genome phylogeny and multi-alignments. Parsnp leverages contextual information provided by multi-alignments surrounding SNP sites for filtration/cleaning, in addition to existing tools for recombination detection/filtration and phylogenic reconstruction.

The tldr is that Parsnp is meant for aligning large sets of very closely related microbial sequences. For example, others in my lab have used parsnp to align over a million SARS-CoV-2 genomes. The only potential concern is that the Parsnp workflow also finishes with a non-optional phylogeny estimation step that takes quite long on such large datasets. Someone should really make it so that you can opt out of that step...

Spoiler/disclaimer, as the current maintainer of Parsnp, that someone is probably me :grimace:

The output of Parsnp is in XMFA. You're probably looking for a fasta alignment, so you can use HarvestTools to convert them:

```
harvesttools -x sars_align/parsnp.xmfa -M sars_align/parsnp.maf

```

Both parsnp and harvest-tools are available on bioconda, and I'd highly recommend installing them that way.

---
Tags: protein-structure, multiple-sequence-alignment
---


# Question
Title: Independent Subset of Rectangles Problem for Nonoverlapping local alignments
I'd like to implement a graph-based algorithm for DNA comparison, more especially by solving a problem that can be formulated as a Maximum Weighted Independent Set problem.

I have found the article in reference \[1\], especially part 2.2.2 "The Independent Subset of Rectangles (IR) Problem", but it is not very explicit on the way to build the graph G given 2 DNA sequence

So I'm looking for any indication, article, or description on the way to construct the conflict graph cited in article reference \[1\] and \[2\].

I'm not from the biology world but the computer world.

Thank you very much for your help.

\[1\] Efficient combinatorial algorithms for DNA sequence processing: https://www.cs.uic.edu/~dasgupta/resume/publ/papers/book-chapter-dasgupta-kao.pdf

\[2\] Nonoverlapping local alignments (weighted independent sets of axis-parallel rectangles) : https://core.ac.uk/download/pdf/82057963.pdf

# Answer
> 1 votes
Consider the below image from a recent review on whole genome alignment:

The lines in the figure are local alignments (sometimes referred to as *anchors* in literature). These lines can be maintained a number of ways, but often via pairwise genome alignment methods like LastZ, MUMmer, and Last (there are many many more).

Often times, one wants to pull out the "best set" of local alignments. In some cases, algorithms will use the "single copy" heuristic, which constrains that each nucleotide in the query can be aligned to at most one nucleotide from the reference. In this case, you can create a bounding rectangle along each of the lines, and two rectangles are independent if there is no point within either rectangle that is on the same $x$ or $y$ coordinate as any point in the other rectangle. You can use this to set your constraint graph.

However, it is worth noting that two rectangles may overlap just a tiny bit. In this case, it would be ideal to chop the overlapping part off of either rectangle, that way you can select both and not just one. This requires a bit more methodology, but a crude way would be to just look for such cases and for each one, chop the overlap off of the less scoring rectangle.

---
Tags: phylogenetics
---


# Question
Title: Evolutionary speaking, what are the constraints on having new mutations?
Evolutionary speaking, some mutations lead to better fitness of an organism and its adaptation to the environment changes, but I was wondering if there are some constraints on having new mutations. For example, if a gene is needed for survival, I expect to see it more resistant to having new mutations. So, my question is what are these constraints, especially for bacteria?

# Answer
This is an undergraduate/graduate lecture type question and the answer varies: traditional Darwinists believed a reasonable proportion of mutations would result in beneficial adaptation, neutral theory however would state the vast majority of mutations have no effect on phenotype and only an absolute minority result in beneficial adaptation. It depends on which bacteria are the focus of investigation and the selective constraint (I would assume its antibiotics).

Once the bacterial species and selective constraint are described its a straight-forward empirical investigation using molecular mathematics because the mutations driving antibiotic resistance are all known.
> 0 votes

---
Tags: genome, bacteria, positive-selection
---


# Question
Title: use same output in two processes in nextflow dsl2
This is my workflow:

```
pre_align()
pre_align.out.single_fastqs.view()
get_fq_info(pre_align.out.single_fastqs)

align_bwa(get_fq_info.out.fq_info)

align_bowtie2(get_fq_info.out.fq_info)

```

where I want to use the same output from `get_fq_info` as input for either `align_bwa` and align\_bowtie2\`.

As it's shown here **only align\_bwa is executed**. If I comment `align_bwa` THEN `align_bowtie2` gets executed.

How can I fix this?

Thank you!

# Answer
> 2 votes
I have solved this problem in the past by duplicating the output channel of the desired process in order to feed the same output into multiple downstream processes.

# Answer
> 1 votes
As far as I know, DSL2 should be able to handle this. Have you specified an output qualifier for your process `get_fq_info` (e.g. `path`)? I ran into a similar issue when I forgot to specify it.

For example, this will result in the error `process align_bowtie2 declares 1 input channel but 0 were specified`:

```
process get_fq_info {
   input:
   path reads

   output:
   "output.txt"

   shell:
   """
   myprocess.sh ${reads} > output.txt
   """
}

```

And is probably solved by this:

process get_fq_info {
   input:
   path reads

   output:
   path "output.txt"

   shell:
   """
   myprocess.sh ${reads} > output.txt
   """

```

---
Tags: nextflow
---


title: VCF file inspector app with GUI

text:
# Question
Title: VCF file inspector app with GUI
Does anyone know a good, free-to-use VCF file inspector app (Mac OS, with GUI) where I can browse the content of a VCF file similarly to bamseek\*?

I would like to browse SNP annotation split by columns, possibly even filter it.

'\* Bamseek does not work for me, fails at open file. It anyways seems to be an abandonware.

Example table view from an expensive commercial software:

# Answer
> 1 votes
I wrote vcf2table: http://lindenb.github.io/jvarkit/VcfToTable.html

```
bcftools view in.bcf |  java -jar dist/vcf2table.jar | more

```

and http://lindenb.github.io/jvarkit/SwingVcfView.html

```
java -jar dist/swingvcfview.jar indexed.vcf.gz

```

# Answer
> 2 votes
OK it's a command line program, not a GUI, but `bcftools query` and `bcftools view` can easily be used to to view elements of a vcf/bcf file and are (probably) a lot more flexible and scalable than a GUI. Examples from here:

```
$ bcftools query -f '%CHROM %POS %REF %ALT\n' file.bcf | head -3
1 13380 C G
1 16071 G A
1 16141 C T

```

```
$ bcftools query -f '%CHROM %POS  GTs:[ %GT]\t PLs:[ %PL]\n' file.bcf | head -3
1 10234  GTs: 1/1 1/1    PLs: 28,3,0 29,3,0
1 10291  GTs: ./. 1/1    PLs: 0,0,0 28,3,0
1 14907  GTs: 0/1 0/1    PLs: 8,0,17 26,0,48

```

---
Tags: vcf, snp
---


title: Counting reads for each biotype

text:
# Question
Title: Counting reads for each biotype
I have bam files. I want to find the total read counts associated with all the biotypes, eg snRNA,rRNA,tRNA mRNA,scRNA,snoRNA etc. I can use ht-seq count to get read counts for the genes, but is there a tool which can directly sum up counts for each of the above category. I have a gencode hg38 gtf file as my annotation file.

Thanks

# Answer
Maybe there's one tool that does everything, but here's one approach that uses a mix of Unix tools and BEDOPS.

At the start, one thing that may be necessary is to fix GTF files that lack an attribute required by specification:

```
$ awk '{ if ($0 ~ "transcript_id") print $0; else print $0" transcript_id \"\";"; }' annotations.gtf > annotations.fixed.gtf

```

This may not be an issue for your GTF file.

With a correctly-formatted GTF, we can filter it by `gene_biotype` attribute, e.g., `snoRNA`, and convert those filtered annotations to a sorted BED file with `gtf2bed`:

```
$ NEEDLE="snoRNA"
$ grep -v '^#' annotations.fixed.gtf | awk -vneedle=${NEEDLE} '{ match($0, /gene_biotype "([a-zA-Z_]+)"/, a); if ((needle==a[1]) && ($8=="gene")) { print $0; } }' | gtf2bed - > ${NEEDLE}.bed

```

Secondly, we can convert the BAM files to sorted BED files via `bam2bed`:

```
$ bam2bed < reads.bam > reads.bed

```

Once we have `snoRNA.bed` and `reads.bed`, we can use `bedmap --count` to count the number of reads which overlap each snoRNA-annotated gene from the original annotations:

```
$ bedmap --count snoRNA.bed reads.bed > answer.txt

```

One complication is that the reads should have a chromosome name scheme identical to that of your GTF file, so that mapping can be done between elements on the same chromosome. Mixing Ensembl and UCSC sourced files, for example, could require an extra step to add/remove `chr` prefix.

In any case, each line of `answer.txt` is the number of reads that overlap a snoRNA-typed annotation, not the overall sum.

To sum all values, you can add to the `bedmap` statement an `awk` command that runs a simple accumulator:

```
$ bedmap --count snoRNA.bed reads.bed | awk 'START{s=0;}{s+=$0;}END{print s;}' > sum.txt

```

The file `sum.txt` will have the sum of reads overlapping annotations.

This procedure could be repeated for other gene biotype categories, `snRNA`, `rRNA`, etc. by making a script that swaps in a category name from a list of such, rerunning the GTF conversion and filtering step, and re-running the `bedmap` step on the new annotation subset.
> 1 votes

# Answer
You may be interested in a Python library I developed called PyReference

```
# Install
python3 -m pip install pyreference
# Process GTF into a JSON file
pyreference_gff_to_json.py --gtf gencode.hg38.gtf
# Create a setup file (see website for details)
# Run against your BAM file
pyreference_biotype.py ${BAM_FILE}

```

This generates a graph:

As well as a CSV of read counts per biotype for each length. You can sum up the read lengths for a total.
> 0 votes

---
Tags: sam, reads, gtf
---


# Question
Title: How fo force nextflow to repeat a process until all values in a particular channels are used up BUT a single value from another channel is needed
I am trying to find a hack for getting a process to run until all emissions in a channel are used up. The problem is that the process also takes as input a second channel that only has a single emission. Therefore, the process only runs once, but I need it to repeatedly use the single emission channel for all the emissions in the other channel. I imagine this is not too uncommon a problem, but I have not found a solution.

To be more clear, here is the actual bioinformatics I am trying to do: I have split fastqs into files for multiple samples and emit these files into a channel. I copy this channel because I want to use the files again later. I then build an assembly from the combined files which I index with bwa. I next pass the indexed assembly (index channel) into a new process that also takes as input ALL of the original fastq files (the ones I put into a duplicate channel earlier). The goal is to align each of the fastq files against the genome. However, there are only a few emissions in the index channel so the process doesn't repeat for all the items in my fastq channel.

To aid you all here is a simple example that may help, but take note that my situation doesn't use value channels so is a bit trickier to deal with.

```
exampleChannel_multiplethings = Channel.from( [1, 3, 5, 7, 9] )
exampleChannel_onething = Channel.from( [1] )

process foo {
  input:
  val x from exampleChannel_multiplethings
  val y from exampleChannel_onething

  output:
  file 'x.txt' into foo_out

  """
  echo $y $x > x.txt
  """
}

```

This outputs "1 1" into the x.txt file. This is expected behavior, but I need to circumvent it somehow so that the process is repeated for all the items in the exampleChannel\_multiplethings channel. Important note, if I was dealing with single values I could probably use one of the combine operators, but in reality I am trying to pass a handful of files for an indexed genome into the process. I haven't found a combine operator that works, but hopefully I have missed something!

Thanks for any assistance

Actual code (overly lengthy to paste here) can be found at: https://github.com/JHarrisonEcoEvo/nextflow\_gbs/blob/main/gbs\_nextflow

# Answer
I discovered a solution using a combination of `.collect()` and `each`.

If one collects from a path the result is a value channel of paths. Value channels can be used multiple times, whereas queue channels cannot.

```
ch_collected = Channel.fromPath( "somepath" ).collect()

```

`ch_collected` can be used multiple times but if one tries to input it to a process using the usual file..from syntax one will get all the files input to the process at once. Instead, one must use the `each` option during input. E.g.,

```
  input:
    each fastq_file from ch_collected

```

I hope this helps someone.
> 4 votes

---
Tags: bwa, nextflow
---


# Question
Title: Is there a JSON-based genomic feature format?
There's a variety of formats for storing genomic features (e.g. Genbank, GFF, GTF, Bedfile) but all of the ones that I'm familiar with use either a custom format (Genbank) or a CSV with defined column conventions and space for additional data (GFF, GTF, Bedfile).

Is there a feature format that is JSON based and adheres to a defined schema?

# Answer
> 5 votes
JVCF is a JSON schema for describing genetic variants. BioJSON is a JSON schema for multiple sequence alignments. The output of mygene.info's gene annotation web service is JSON-formatted.

You might build one of your own schemas from any or all of these, as inspiration for structure and content.

Consider that most genomics formats are historically tabular, because they are not hierarchical, because they originated prior to the advent of JavaScript as a common programming language (and subsequent use of JSON to serialize data), and because a lot of stream-oriented UNIX tools read in data and process it in ways that are natural to tabular structure.

Some tabular results can be hierarchical, but use inefficiencies in field values to apply hierarchy. JSON is inefficient for storing tabular data, but it is absolutely wonderful for hierarchical data. This difference is perhaps worth considering.

One could think of storing a gene as a master object containing various haplotypes, exons, introns, features, associated variants, associated TF binding sites, gene ontologies, etc. Ordering or prioritization of such items would be possible with the use of sorted lists within a gene object. The format is as extensible as it is a matter of adding a property, list, or object to the schema.

To gain adoption or "buy-in", build tools which publish it to users, and tools that consume what users get in order to do things, along with conversion and validation utilities.

# Answer
> 1 votes
BioJS has capacity to read and write GFF as streams, it might have something of interest such as this: https://github.com/biojs-io/biojs-io-gff

# Answer
> 1 votes
I've written a Python library called PyReference which converts a GTF / GFF (RefSeq and Ensembl) into a gzipped JSON file.

You can just use this if you like but there's also a Python wrapper around the JSON, which allows you to write genomics code more naturally.

Python isn't known for being fast but the library function for reading a JSON file is highly optimised - the following takes less than 4 seconds on my laptop:

```
import numpy as np
import pyreference

reference = pyreference.Reference()

my_gene_ids = ["MSN", "GATA2", "ZEB1"]
for gene in reference[my_gene_ids]:
    average_length = np.mean([t.length for t in gene.transcripts])
    print("%s average length = %.2f" % (gene, average_length))
    print(gene.iv)
    for transcript in gene.transcripts:
        if transcript.is_coding:
            threep_utr = transcript.get_3putr_sequence()
            print("%s end of 3putr: %s" % (transcript.get_id(), threep_utr[-20:]))

```

---
Tags: features, genomefeatures
---


# Question
Title: How to remedy a DESeq2 collapsing technical replicates error?
**Goal:** To ensure "the sum of the counts for \[my samples\] is the same as the counts in the \[samples\] columns in `ddsColl`" after collapsing replicates using DESeq2.

Before performing differential expression (DE) analysis on my RNA-seq data stored in a DESeqDataSet object, I need to collapse the read counts observed for each of my technical replicates. This should be simple enough using the DESeq2 manual and vignettes available as reference.

This is how you collapse replicates:

```
ddsColl <- collapseReplicates(dds, dds$sample, dds$run)
# examine the colData of the collapsed data
colData(ddsColl)
# check that the sum of the counts for "sample1" is the same
# as the counts in the "sample1" column in ddsColl
matchFirstLevel <- dds$sample == levels(dds$sample)[1]
stopifnot(all(rowSums(counts(dds[,matchFirstLevel])) == counts(ddsColl[,1])))

```

However, when I try this, I get the following error with no apparent resolution reported in the manual or on online forums:

```
Error: all(rowSums(counts(dds[, matchFirstLevel])) == counts(ddsColl[,  .... is not TRUE

```

**How might I circumvent this error and continue on with the analysis?**

---

**Input files:**

```
> head(cts)
            KO1_P1_L001.Counts KO1_P1_L002.Counts KO1_P2_L001.Counts KO1_P2_L002.Counts KO1_P3_L001.Counts KO1_P3_L002.Counts KO2_P1_L001.Counts KO2_P1_L002.Counts KO2_P2_L001.Counts KO2_P2_L002.Counts KO2_P3_L001.Counts KO2_P3_L002.Counts WT_P1_L001.Counts WT_P1_L002.Counts WT_P2_L001.Counts WT_P2_L002.Counts WT_P3_L001.Counts WT_P3_L002.Counts
DDX11L1                      1                  0                  1                  1                  2                  1                  0                  1                  2                  1                  3                  1                 7                 2                 1                 1                 0                 0
WASH7P                     124                144                128                151                102                118                 39                 41                103                 87                105                106               125               120                13                26               104               136
MIR6859-1                    8                  6                  4                  5                  4                  6                  2                  2                  7                  2                  7                  6                 6                 6                 2                 2                 2                 9
MIR1302-2HG                  1                  4                  0                  0                  0                  2                  2                  2                  2                  2                  0                  4                 1                 2                 0                 3                 2                 2
MIR1302-2                    0                  0                  0                  0                  0                  0                  0                  0                  0                  0                  0                  0                 0                 0                 0                 0                 0                 0
FAM138A                      0                  0                  0                  0                  0                  0                  0                  0                  0                  0                  0                  0                 0                 0                 0                 0                 0                 0

> print(coldata)
                   condition sample           run
KO1_P1_L001.Counts    longKO     l1  long_P1_L001
KO1_P1_L002.Counts    longKO     l1  long_P1_L002
KO1_P2_L001.Counts    longKO     l2  long_P2_L001
KO1_P2_L002.Counts    longKO     l2  long_P2_L002
KO1_P3_L001.Counts    longKO     l3  long_P3_L001
KO1_P3_L002.Counts    longKO     l3  long_P3_L002
KO2_P1_L001.Counts   shortKO     s1 short_P1_L001
KO2_P1_L002.Counts   shortKO     s1 short_P1_L002
KO2_P2_L001.Counts   shortKO     s2 short_P2_L001
KO2_P2_L002.Counts   shortKO     s2 short_P2_L002
KO2_P3_L001.Counts   shortKO     s3 short_P3_L001
KO2_P3_L002.Counts   shortKO     s3 short_P3_L002
WT_P1_L001.Counts         WT     w1    wt_P1_L001
WT_P1_L002.Counts         WT     w1    wt_P1_L002
WT_P2_L001.Counts         WT     w2    wt_P2_L001
WT_P2_L002.Counts         WT     w2    wt_P2_L002
WT_P3_L001.Counts         WT     w3    wt_P3_L001
WT_P3_L002.Counts         WT     w3    wt_P3_L002

```

**My workflow:**

```
> dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~ condition)
> ddsColl <- collapseReplicates(dds, dds$sample, dds$run)
> colData(ddsColl)
DataFrame with 9 rows and 4 columns
   condition      sample           run               runsCollapsed
    <factor> <character>   <character>                 <character>
l1    longKO          l1  long_P1_L001   long_P1_L001,long_P1_L002
l2    longKO          l2  long_P2_L001   long_P2_L001,long_P2_L002
l3    longKO          l3  long_P3_L001   long_P3_L001,long_P3_L002
s1   shortKO          s1 short_P1_L001 short_P1_L001,short_P1_L002
s2   shortKO          s2 short_P2_L001 short_P2_L001,short_P2_L002
s3   shortKO          s3 short_P3_L001 short_P3_L001,short_P3_L002
w1        WT          w1    wt_P1_L001       wt_P1_L001,wt_P1_L002
w2        WT          w2    wt_P2_L001       wt_P2_L001,wt_P2_L002
w3        WT          w3    wt_P3_L001       wt_P3_L001,wt_P3_L002

> matchFirstLevel <- dds$sample == levels(dds$sample)[1]
> matchFirstLevel
logical(0)

> stopifnot(all(rowSums(counts(dds[,matchFirstLevel])) == counts(ddsColl[,1])))
Error: all(rowSums(counts(dds[, matchFirstLevel])) == counts(ddsColl[,  .... is not TRUE

```

From the output of `dds` and `ddsColl`, it seems the number of columns decreased from 18 to 9, as expected. The only difference I see between my DESeq2 workflow output is the result of `matchFirstLevel` (see below for expected result). Why might this be for my dataset, but not the example dataset?

**Example workflow:**

```
> dds <- makeExampleDESeqDataSet(m=12)
> dds$sample <- factor(sample(paste0("sample",rep(1:9, c(2,1,1,2,1,1,2,1,1)))))
> dds$run <- paste0("run",1:12)
> ddsColl <- collapseReplicates(dds, dds$sample, dds$run)
> colData(ddsColl)
DataFrame with 9 rows and 4 columns
        condition   sample         run runsCollapsed
         <factor> <factor> <character>   <character>
sample1         B  sample1       run10   run10,run11
sample2         A  sample2        run5          run5
sample3         A  sample3        run6          run6
sample4         A  sample4        run1     run1,run4
sample5         B  sample5        run7          run7
sample6         B  sample6        run9          run9
sample7         A  sample7        run3    run3,run12
sample8         B  sample8        run8          run8
sample9         A  sample9        run2          run2

> matchFirstLevel <- dds$sample == levels(dds$sample)[1]
> matchFirstLevel
 [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE FALSE

> stopifnot(all(rowSums(counts(dds[,matchFirstLevel])) == counts(ddsColl[,1])))

```

# Answer
`sample` is a character vector, so `levels(dds$sample)` is returning `null`. Presumably you mean to use `dds$sample[1]` here instead, which will likely solve the problem.
> 1 votes

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


# Question
Title: Mounting a folder from the host to a Docker, while running it as Cromwell runtime
I have a pipeline, which is written in `WDL` & runs with the help of `Cromwell`. I use `Docker` for containerization.
However, my container is dependent on huge (\>100GB databases), which cannot be shipped with the image.
Is there an option to mount a directory from the host filesystem to a `Docker` while running `Cromwell` workflow? I found this functionality in `snakemake`, therefore I assume it is common, but I haven't found it in the docs yet.

# Answer
> 1 votes
Not an expert of WDL but looking at Cromwell documentation, there is a block in the configuration of `Cromwell` called `submit-docker`:

> you'll need to modify the config block inside backend.providers in your Cromwell configuration. In particular, this block contains a key called submit-docker

In this example file, on line `420`, it is possible to specify the mount option which is `-v /path/you/need/to/mount:/path/in/docker/container` (usually they are the same) to Docker command line. Setting it like that should work.

---
Tags: wdl, workflow-management, docker
---


# Question
Title: How to format my file using linux tools?
I have an imputed dosage for chr 22 file. its content is:

```
#[1]CHROM       [2]POS  [3]REF  [4]ALT  [5]HG00096_HG00096      [6]HG00097_HG00097      [7]HG00099_HG00099      [8]HG00100_HG00100
22      16051249        T       C       0.0     1.0     1.0     0.0    0.0     1.0     0.0     0.0     1.0     0.0

```

I want to remove `[]` and remove anything coming after `_` from this file as well as instead of separate `22` `16051249` `T` and `C` i want to write:

```
22_16051249_T_C_b37

```

So the final file should look like:

```
Id      HG00096 HG00097 HG00099 HG00100
22_16051249_T_C_b37     0       1       1       0       1       0       0       1       0

```

Can anyone please help me out with this.

# Answer
> 0 votes
Assuming that the third line is the continuation of the second line. Here's a one (long) liner:

```
cat <( \
      paste -d '\t' \
         <(echo "Id") \
         <(head -1 test.tsv \
      | cut -f5- \
      | sed -e '$s/\[[[:digit:]]\+\]//g; s/_HG[[:digit:]]\+//g') ) \
    <( \
      paste -d '\t' \
        <(awk 'NR>1 {print $1"_"$2"_"$3"_"$4"_b37"}' test.tsv) \
        <(awk 'NR>1 {print $0}' test.tsv | cut -f5- | sed 's/\.[[:digit:]]\+//g')) \
> output.tsv

```

---
Tags: file-formats
---


title: imputed dosage values for vcf files

text:
# Question
Title: imputed dosage values for vcf files
I have one more question related to genotype file. I submitted the job in michigan imputation server.

and my imputed dosage for chr 22 file is:

> 22 16050435 22:16050435:T:C T C . PASS AF=0.00098;MAF=0.00098;R2=0.23769;IMPUTED GT:DS:HDS:GP 0|0:0.002:0.001,0.000:0.998,0.002,0.000 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0.001:0.000,0.001:0.999,0.001,0.000 0|0:0:0,0:1,0,0 0|0:0.001:0.001,0.000:0.999,0.001,0.000 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0.001:0.001,0.000:0.999,0.001,0.000 0|0:0.002:0.000,0.002:0.998,0.002,0.000 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0.001:0.000,0.001:0.999,0.001,0.000 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0.002:0.000,0.002:0.998,0.002,0.000 0|0:0:0,0:1,0,0 0|0:0.001:0.001,0.000:0.999,0.001,0.000 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0

Its coming like this.

wherein the test data provided in your predixcan tutorial looks like this:

> 22\_16494187\_C\_A\_b37 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0

Also it states that the imputed dosage should be encoded on a 0-2 scale representing the number imputed or number of the effect alleles the sample posseses. Do anyone know should I convert my vcf files to the desired input genotype file to train the model. I don't understand in my vcf files which column is representing dosage values Thank you.

# Answer
When i was still working in the field of genotype imputation, a dosage was the expectation (average) of the (posterior) genotype probabilities, coded as GP in the OP's VCF. I once contributed code to bcftools to convert genotype probabilities stored in the GP field to dosages. It looks like that code may now live in a bcftools plugin called "dosage": https://samtools.github.io/bcftools/howtos/plugins.html
> 1 votes

---
Tags: vcf, snp, genotyping, imputation, chromosomes
---


# Question
Title: What is the cutoff distance in this algorithm?
In the following article:

What are the minimum and maximum cutoff distances in this algorithm?

Also, what does it mean by:

>

I.e. where does (2x43)x3 come from?

# Answer
The 43 features are listed in the supplementary materials on page 4. For example,

```
Feature Type Atoms Involved
1 Distance i, i-5
2 Distance i, i-4
3 Distance i, i-3
4 Distance i, i-2
...
11 Angle i-1, i, i+1
12 Dihedral i, i+1, i+2, i+3

```

Are some of them; so Feature 11 is the angle between atom(i-1), atom(i), and atom(i+1). Since it says these were calculated for both the C-alphas and pseudo-centers, this gives 2 * 43. Then this is multiplied by three as mentioned in your quoted paragraph.

These features (measurements, almost) are numbers, not booleans - so I do not think that they use distance cutoffs as such, although the random forest method is presumably learning various cutoffs and encodes those in the structure of the forest.
> 2 votes

---
Tags: proteins, algorithms
---


# Question
Title: How does DESeq2 "collapseReplicates()" function work on read counts data?
Comparing read counts from an RNA-seq experiment for two select genes *before* and *after* using DESeq2's `collapseReplicates()` and `plotCounts()` functions yields interesting results:

**Before `collapseReplicates()` and `plotCounts()`:**

| Geneid | foo1.1 | foo1.2 | foo2.1 | foo2.2 | bar1.1 | bar1.2 | bar2.1 | bar2.2 | baz1.1 | baz1.2 | baz2.1 | baz2.2 | baz3.1 | baz3.2 |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| WASH7P | 6 | 5 | 0 | 2 | 1 | 1 | 8 | 5 | 0 | 0 | 0 | 0 | 0 | 0 |
| SOX3 | 1880 | 1861 | 1950 | 2055 | 1189 | 1181 | 2415 | 2482 | 3887 | 3810 | 1851 | 1738 | 3217 | 3406 |

**After `collapseReplicates()` and `plotCounts()`:**

| Geneid | foo1 | foo2 | bar1 | bar2 | baz1 | baz2 | baz3 |
| --- | --- | --- | --- | --- | --- | --- | --- |
| WASH7P | 9.877191 | 2.279384 | 3.478891 | 11.613875 | 0.500000 | 0.500000 | 0.500000 |
| SOX3 | 3189.598 | 3563.717 | 3530.486 | 4187.011 | 6473.122 | 7991.460 | 5390.721 |

---

**Note:** In the above comparison, there are two (2) technical replicates (1.1 + 1.2, 2.1 + 2.2, etc.) for each biological replicate (foo1, foo2, etc.) for each of three (3) conditions (foo, bar, and baz).

---

Comparing the tables above, it appears as though DESeq2 is *NOT* taking the average *or* sum of columns being collapsed.

It is also curious - and mildly concerning - that some very low expression (i.e., 0.5 counts) is reported for genes in the matrix of collapsed replicates when, in the original count matrix, zero (0) reads were counted as 'mapped to that gene'.

So, how does `collapseReplicates()` "combine counts into single columns of the count matrix" as is described in the DESeq2 vignette?

---

Here is the code to collapse replicates and retrieve the number of read counts for a specific gene (e.g., *WASH7P*) in a dataframe to be used in a count plot, getting the "after" results shown above:

```
dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~ condition)
ddsColl <- collapseReplicates(dds, dds$sample, dds$run)
keep <- rowSums(counts(ddsColl)) >= 10
dds <- ddsColl[keep,]
dds$condition <- relevel(dds$condition, ref = "baz")
dds <- DESeq(dds)
countsdf <- plotCounts(dds, gene="WASH7P", intgroup="condition", returnData=TRUE)

```

# Answer
> 1 votes
The DESeq2 function `collapseReplicates` sums the counts for the technical replicates.
Here is the code reference: github.com/mikelove/DESeq2/blob/master/R/helper.R#L186

OPs actual confusion was with the DESeq2 function `plotCounts` which by default normalizes count data and adds a pseudocount of 0.5 for plotting on log2 scale.

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


title: PDBQT file format

text:
# Question
Title: PDBQT file format
I am trying to understand the source code of the AutoDock Vina, but as I am new to the field of Bioinformatics it is a little hard to understand. It seems like they are using a tree data structure to represent the ligand and the receptor. I cannot understand what nodes and branches are in this data structure (Heterotree). I think in order to understand it I should first understand what information the PDBQT file contains. Can you tell me what branch and root is in PDBQT file format?

# Answer
> what information the PDBQT file contains

The PDBQT format is described here: https://autodock.scripps.edu/wp-content/uploads/sites/56/2021/10/AutoDock4.2.6\_UserGuide.pdf —ignore the expired SSL certificate warning as I am pretty sure the Scripps Institute is not a criminal organisation.

> what branch and root is in PDBQT file format

It is rather self evident but with a twist and as you have guessed it has to do with trees in graph networks. Lines of text \<-\> graph networks describing a molecule is a common nuisance in compchemistry. Normally the aim is connectivity, but here it is torisons. But the former is a nice example, a small molecules can be represented with SMILES (strings), where say isopropanol is `CC(O)C`, in which you have a branch going off in round brackets, whereas everything else is assumed connected sequentially. In the case of the PDBQT format you are describing what atom subsets of the molecule rotate together and actually not how they connect as given away by the fact that you don't have a special syntax for ring closures (in a SMILES you write cyclopropane as `C1CC1`, where the third carbon connects to C1, which was declared first. This is because the **topology** is dealt with elsewhere.

### Resources

AutoDock 4 and AutoDock Vina are divergent pieces of software, but they have several common elements like the PDBQT format, so you get double the documentation (and papers to read)!
> 1 votes

---
Tags: proteins, pdb, docking
---


title: ClustalOmega alignment

text:
# Question
Title: ClustalOmega alignment
I am performing multiple sequence alignment using Clustal Omega for 600 RNA covid-sequences, i.e. cDNA on Genbank, with input characters around 30574 characters for each sequence.

I am running it on windows cmd. I defined parameters in input of maxseqle = 37000, however, the output each time gives duplicated length around 75000 charcters per each sequence.

```
clustalo.exe -i Allseq.fasta --is-profile --use-kimura --seqtype DNA --maxseqlen 37000 --threads 8 -o myclustalv3.fasta

```

How can I solve this problem, what will I define in input parameters. This problem occurs with large no of sequences 600 sequences. In test of 40 seq, it gave me normal output length around 32000 characters but when align all 600 sequences together, it gave this duplicated length results.

# Answer
Have you had a look at your alignment? I personally use Jalview but any alignment viewer should work. Colour \> Percentage Identity (or Nucleotide) should help with visualization.

Does the alignment look reasonable (and gaps are just inserted at the end for some reason), or is it pretty messy? If it looks messy, is it just a few outlier sequences that cause the problem, or does it look like you have two main groups of sequences that align between themselves but not with the other? In the latter case, could it be possible that some of your initial sequences were inputed in the reverse orientation? I'm not sure about the command line implementation of Clustal Omega, but I know that at least the versions I use don't automatically reverse-complement sequences to be aligned.
> 2 votes

# Answer
Firstly its good you've done some proofing of your output wherein you spotted,

> 75000

The alignment output is approximately double the genome size, just over double the input.

An off-question point, I assume the computer you are working on has at least 8-cores because thats what you are requesting.

Clustal omega is a very solid algorhithm and has been around for a long time in its current omega version. The error you are encountering is an input sequence bug. @Laura's suggest is correct and is formal practice: you must manually verify bioinformatics output and in this case an alignment viewer is the correct proceedure.

Thereafter I suspect you would find the answer, such as omitting carriage returns in the fasta format, seeing taxa names caught up post 30K.

If you run,

```
wc ./myclustalv3.fasta # output file
wc ./Allseq.fasta  # input file
wc ./testgenome.fasta

```

You will probably find that the number lines per file does not agree between the input and output format OR the number of lines of both input and output is \<\<600 or \<\<1200 (600 x2 ), i.e. the expected number of sequences. This points to parsing errors within the input data. Fasta has a line for the sequence identifier and a line for the sequence. It is using \> for the sequence identifier, but will use a carriage return for sequence.

Again manually assessing an alignment is baseline proceedure regardless of output.

An example (guess) about what might have happened is that every other input has omitted "\>", omega will then see the first "\>" then read everything until it sees another "\>" as sequence data. Thus its doubled the sequence length and dragged a sequence identifier exactly in the centre of the doubled up genomes.

---

An issue with fasta is there is no internal proofing of the format. In other formats the number of sequences and their length is specified within the file, so it is easy to internally verfiy the data prior a calculation. Code-wise it is known as 'setting' and you can't do that with fasta.
> 0 votes

---
Tags: sequence-alignment, multiple-sequence-alignment
---


# Question
Title: Getting hierarchy of cell populations with Drop-seq data
I plan on collecting some drop-seq data from brain tissue (from mice and humans). Using only the 75-85% of genes that have 1:1 orthology between mice and humans, I'd like to cluster the cells by transcriptomic similarity and get a dendrogram to use for cell type classification. There are tons of tools for processing and clustering single cell RNA-seq data (https://www.scrna-tools.org/tools), but I get the impression that most of them were not designed with Drop-seq in mind. Anyone have any insight into the pros and cons of the best tools for Drop-seq analysis?

# Answer
This article uses the freely available R package dropbead for filtering and then Seurat to perform a principal component analysis that groups together affine transcriptomes. It could be what you are looking for.
> 4 votes

# Answer
SC3 is part of the Bioconductor library for R. It is based on PCA and spectral dimensionality reductions and utilizes k-means. In the introductory paper the authors perform analysis on a Drop-seq dataset of more than 44,000 cells. SC3 outputs several useful plots, including a consensus matrix, expression panel, and a marker genes heatmap. All of these utilize the hierarchical clustering.

I can't speak from personal experience but the package seems very user friendly. Being a Bioconductor package means it has lots of support and excellent documentation. I am getting ready to use it myself for a small Drop-seq dataset of olfactory epithelial cells.

---

SC3: consensus clustering of single-cell RNA-seq data : Nature Methods : Nature Research. Available at: https://www.nature.com/nmeth/journal/v14/n5/full/nmeth.4236.html. (Accessed: 31st October 2017)
> 1 votes

# Answer
You can try a hierarchical clustering tool HGC (hierarchical graph clustering).

You may get cell population at different hierarchy by choosing multiple cutting heights. See the paper for more details.
> 1 votes

# Answer
The current version of Seurat and Scanpy can both handle Drop-seq data. If you want to include more low-reads cells you can loosen the filtering thresholds in either of these packages. These two packages really can handle both low-read-high-cell data and high-read-low-cell data.

This ENCODE paper used Seurat to process both C1 Fluidigm data (\>4000 genes per cell) and 10X data (\>1000 genes per cell only).

This human limb paper used Seurat to process Drop-seq data, although they didn't annotate the cell types well. I used Seurat to reprocess their data and it works much better than what the paper presents.

I guess the key is parameter tuning.
> 0 votes

---
Tags: rna-seq, clustering, single-cell
---


# Question
Title: #1 Parsing trouble java.lang.NumberFormatException in GSEA 4.2.1
I am getting the following error message when I tried to load .TXT data file to GSEA 4.2.1. I tried different way to change the file name as .GCT but not working. Any suggestion? Many thanks.

Full Error Message:

```
There were errors: ERROR(S) #:1
Parsing trouble
java.lang.NumberFormatException: ...

---- Stack Trace ----
# of exceptions: 1
------For input string: "description"------
java.lang.NumberFormatException: For input string: "description"
    at java.base/jdk.internal.math.FloatingDecimal.readJavaFormatString(Unknown Source)
    at java.base/jdk.internal.math.FloatingDecimal.parseFloat(Unknown Source)
    at java.base/java.lang.Float.parseFloat(Unknown Source)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.AbstractParser.parseStringToFloat(AbstractParser.java:266)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.AbstractParser.parseFieldsIntoFloatArray(AbstractParser.java:371)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.AbstractParser.parseTextMatrixToDataset(AbstractParser.java:293)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.TxtDatasetParser.parse(TxtDatasetParser.java:109)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserFactory.readDatasetTXT(ParserFactory.java:201)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserFactory.read(ParserFactory.java:760)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserFactory.read(ParserFactory.java:737)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserWorker.doInBackground(ParserWorker.java:52)
    at java.desktop/javax.swing.SwingWorker$1.call(Unknown Source)
    at java.base/java.util.concurrent.FutureTask.run(Unknown Source)
    at java.desktop/javax.swing.SwingWorker.run(Unknown Source)
    at java.base/java.util.concurrent.ThreadPoolExecutor.runWorker(Unknown Source)
    at java.base/java.util.concurrent.ThreadPoolExecutor$Worker.run(Unknown Source)
    at java.base/java.lang.Thread.run(Unknown Source)

```

# Answer
I just figured out the error. My .TXT was not correctly formated, following the instruction from their manual under .TXT format worked.
> 2 votes

---
Tags: gsea, parsing
---


# Question
Title: Assigning subcluster idents to original object
I have a scRNA-seq Seurat object I've analyzed, and I noticed that for some of the clusters, there's more than one cell type.

I've created a subset which and run FindClusters again to label the cell more efficiently, and now I want to "paste" the Idents I've assigned in the subcluster to the original object.

How do I do it?

# Answer
> 1 votes
Since you haven't shared any data or code so I will show it using the example given below,

First you need to assign a specific cluster to the cells which are having another type like C\_new.

```
cells_df <- data.frame(row.names=paste0("Cell", 1:10, sep = ""),as.vector(paste0("C",1:10,sep="")))

colnames(cells_df) <- "clusters"

cluster_new <- c("Cell6","Cell2")

cells_df$clusters <- as.character(cells_df$clusters)

for(i in 1:length(cluster_new)) {
j <- cluster_new[i]
cells_df[j,] <- "C_new"
}

head(cells_df)
         clusters
 Cell1       C1
 Cell2    C_new
 Cell3       C3
 Cell4       C4
 Cell5       C5
 Cell6    C_new

```

And if you have already assigned new clusters to the cells then you can skip the above step. Now save the new cluster belongings into the seurat object using the `Idents` function

```
Idents(seurat_obj) <- cells_df$clusters

```

Then you can perform the differential gene expression as follows, for example if you want to find the differentially expressed genes between cluster `C7` vs `C_new` you can do it as follows,

```
FindMarkers(seurat_obj,ident.1 = "C7", ident.2 = "C_new" )

```

# Answer
> 1 votes
Assuming `parentObj` is the original seurat object, and `tmp_subset` is your subset, you can

```
new_metadata <- parentObj@meta.data
tmp_subset_annotation <- tmp_subset@meta.data

# new_metadata
           original_annotation
ATTCGGA-1   TCD4
AAATAAA-1   TCD4-naive
...

# tmp_subset_annotation
           new_cellType
ATTCGGA-1   TCD4-naive
AAATAAA-1   TCD4-naive
...

new_metadata$refined_annotation <- sapply(rownames(new_metadata),
     function(ita) ifelse(ita %in% rownames(tmp_subset_annotation),
            tmp_subset_annotation[match(ita,rownames(tmp_subset_annotation)),"new_cellType"],
            new_metadata[match(ita,rownames(new_metadata)),"original_annotation"]))

```

Basically you iterate each rownames of the original object, and assign the new annotation, if the UMI matches a UMI in the subset, otherwise you keep the original annotation.

---
Tags: scrnaseq, seurat, single-cell, clustering
---


# Question
Title: In a RNA-Seq heatmap should you do Z-score standardisation before clustering the rows/columns or after?
I have made a heatmap using RPKM values from a RNA-Seq dataset using the pheatmap() function in R. I have log2-transformed the data before performing Z-score standardisation of the data. I have also clustered the rows and columns of the heatmap. I have the following code:

```
heatmap_data_log2 %>% pheatmap(color = colorRampPalette(c("blue2","white","red"))(100), scale = "column", cluster_cols = T, clustering_method = "ward.D2", angle_col = 45, fontsize_row = 7.5, fontsize_col = 8, border_color = NA, cutree_rows = 4)

```

I have seen in a video tutorial that the pheatmap() function does the Z-score scaling before doing the clustering of the rows and columns. However, in this online article, it says that *"The Z-scores are computed after the clustering, so that it only affects the graphical aesthetics and the color visualization is improved"*.

As these two sources are giving contradictory information, I was wondering which is better. Should you do Z-score scaling of the gene expression values before doing the clustering or after? Any advice is appreciated.

# Answer
Looking at the source for `pheatmap`, there is a function called `scale_mat` that is used to preprocess and normalize the input matrix, depending on the value of `scale`, which specifies one of either none, row, or column normalization options. Separate and similar functions are used for color value scaling downstream.

I have not watched the video or read the linked article, so I'm not sure this is contradictory, but it may be just small confusion about how normalization is applied (and to what).

Normalization is usually done before clustering to focus cluster construction on signal rows or columns that are not some multiple or scale factor of another — and in referring back to the source code, this is what `pheatmap` does, as well.
> 3 votes

---
Tags: r, heatmap
---


# Question
Title: How do I set my path for the SRA toolkit under macOS 12.1?
I would like to add the SRA toolkit to my path under macOS 12.1.

After downloading the SRA directory, decompressing, and adding to the path, the directory does not seem to be available e.g. I'm attempting to add it to my path with

```
export path=(/Users/ray/Documents/School/bioinformatics/sratoolkit.2.11.2-mac64/bin)

```

and it does not work. For example if I type

```
 `ray@dynamic-xx-xxx-xxx-xx ~ % which prefetch

```

prefetch not found\` macOS gives me prefetch not found. The result of my path is

```
`ray@dynamic-10-125-206-18 ~ % echo $PATH
/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin`

```

Is there a way to add it to my path so it works?

# Answer
Firstly, you are trying to set a path into the Documents directory? Thats like for documents.

---

Anyway there could be a permission issue. The quickest way is simply set a path within .bashrc or .profile whatever into ~/bin and then use

```
ln -si /Users/ray/Documents/School/bioinformatics/sratoolkit.2.11.2-mac64/bin/prefetch ~/bin

```

... for each compiled program.

```
$ prefetch

```

should then work. This creates a symbolic link in `~/bin` so when you type the program in your commandline prompt it will fire, without wiring it via a custom bin. If there is a permission issue a permission error will be quickly flagged.

---

The best way is to use Anaconda because its within the Anaconda collection and simplifies everything on three counts, viz:

* in answer to your question all configurations are automated
* it is easy to automate all updates
* it offers a virtual environment where code can't mess with your system
* whenever you want to download other packages its the same approach

Thus all you have to do is download and install Anaconda or miniconda. From this point everything is done using conda, either:

```
conda install -c bioconda sra-tools
conda install -c bioconda/label/cf201901 sra-tools

```

If one fails the other one will work. Then simply,

```
conda activate

```

SRA is now fully configured and ready to go. To get out of the conda virtual environment,

```
conda deactivate

```

If you use Anaconda alot it is better to make a specific environment for SRA and not use the core environment, but thats another question.

Anaconda installs in two different locations depending on whether you install the GUI version, which makes maintenance and environments a bit easier. In the later case I think its `/opt` (but forget), otherwise the commandline version is `~/anadonda3`.

---

I see your using a virtual environment via `homebrew`. Anaconda is separate, its the same idea, but its better.
> 0 votes

---
Tags: sratoolkit
---


# Question
Title: what are the nodes and edges values in gene regulatory networks?
I am trying to find out how one can using gene expression data can infer gene regulatory network applying graph theory concepts. But I could not find a proper reference that

1. explain how one can get the adjacency matrix for a gene regulatory network? like what are nodes what are the edges? which kinds of parameters one should infer?
2. how one can use gene expression data for this modeling?

Any reference or help regarding these questions is highly appreciated if even it is going to be a keyword.

# Answer
This is an incomplete answer since it's not my field of expertise, but it sounds like Weighted correlation network analysis is what you are referring to.

There is a paper which describes how it can be applied to gene expression data.

Langfelder, P., Horvath, S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9, 559 (2008). https://doi.org/10.1186/1471-2105-9-559
> 3 votes

# Answer
There are many different kinds of "gene regulatory networks", so how to model them with gene expression data depends on what you're trying to model and what type of regulatory mechanisms you want to study.

Pick a random gene in the human genome and use a resource database like StringDB. Let's take the TET1 protein for example. Some edges in that graph represent co-expression that you'd be able to find from protein expression data (e.g. NANOG, DNMT3B, and more are linked to TET1 this way). This type of graph tries to model how different proteins interact with each other inside a cell.

Alternatively, you can look at 3D chromatin interactions from Hi-C data, for example. Nodes in that graph could be regulatory elements like enhancers and promoters, and the edges could be how frequently these elements interact with each other. This type of graph tries to model how DNA elements interact with each other to control expression of nearby genes. This is a different purpose to the one above, so it needs to be modelled and interpreted differently.

Each of these are their own type of "gene regulatory network" which require their own methods for constructing adjacency matrices. For example, you could construct an adjacency matrix by looking at the correlation of expression for each pair of proteins (or transcripts, or genes) across a large set of samples. Or in the promoter-enhancer interaction case, a Hi-C contact matrix is itself a type of adjacency matrix.

Different graph models (and thus different adjacency matrices) will be useful for different purposes. So it might be worth trying to refine what you mean by a "gene regulatory network", what aspect of biology you're trying to study using this network, and how different types of gene expression data may be useful to construct these networks.

John Quackenbush's and Andrea Califano's labs have produced a lot of papers modelling different types of gene regulatory networks, like The Network Zoo and ARACNe. These mostly use gene expression data from RNA-seq measurements, but you can read up on some of these to get some ideas on how these types of models may work.
> 1 votes

---
Tags: genomics, gene-expression
---


# Question
Title: Find the genomic locations of exons of certain human gene
So I know Consensus CDS has the genomic locations of exons, however, only of the exons in the CDS. Everything before or after is not included.

NIH variation viewer has those (kinda-) radio buttons next to the gene name, and when you hover over them, you do get to see the genomic locations of the exons. But this method isn't very scalable.

Does anyone know an easy, scalable way of retrieving the genomic locations of exon for a certain gene?

# Answer
> 2 votes
Perhaps you could use this (assuming `hg38`) as a lookup table:

```
$ wget -qO- https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_39/gencode.v39.annotation.gff3.gz \
    | gunzip -c \
    | grep -wF exon \
    | convert2bed --input=gff --attribute-key="gene_name" - \
    > gencode.v39.exons.bed

```

If you prefer to search for an Ensembl accession ID, you can remove the `--attribute-key="gene_name"` option. Otherwise, this puts the HGNC symbol into the ID column.

The preceding step takes about 30 seconds. Once you have made this lookup file, you can search it for exons of any gene of interest. Using `MYC` as an example:

```
$ awk -v FS="\t" -v OFS="\t" '($5=="MYC")' gencode.v39.exons.bed > answer.bed

```

You can then further filter that by transcript or other attributes, etc.

# Answer
> 1 votes
I'd use biomart on ensembl.org

---
Tags: genome, human-genome, sequence-annotation, exon
---


# Question
Title: R package development: How does one automatically install Bioconductor packages upon package installation?
I have an R package on github which uses multiple Bioconductor dependencies, 'myPackage'

If I include CRAN packages in the DESCRIPTION via `Depends:`, the packages will automatically install upon installation via devtools, i.e.

```
devtools::install_github('repoName/myPackage')

```

This is discussed in Section 1.1.3 Package Dependencies, in Writing R Extensions

Is there a way to streamline this such that packages from Bioconductor are automatically installed as well?

Normally, users install Bioconductor packages via BiocLite, e.g.

```
source("http://www.bioconductor.org/biocLite.R")
biocLite("edgeR")

```

# Answer
> 10 votes
There's a trick to this where one needs to add `biocViews:` to the package Description. That's the only solution I've ever seen to allowing automatic installation of bioconductor dependencies. If you need a couple of examples, then click through the link I posted and scroll down to pull requests referencing that issue, they will generally include the actual example.

# Answer
> 18 votes
As suggested, here’s an example showing the relevant lines from a `DESCRIPTION` file from a CRAN/GitHub hosted project that has Bioconductor dependencies (truncated):

```
Depends:
    R (>= 3.3.0)
biocViews:
Imports:
    methods,
    snpStats,
    dplyr

```

The relevant bit is the empty `biocViews:` declaration, which allows the Bioconductor dependency {snpStats} to be automatically installed.

# Answer
> 2 votes
I have a counter example to other answers. BiocViews is not needed to install Bioconductor dependencies on CRAN. I have a package on CRAN that suggests Bioconductor packages and they are found without using BiocViews.

If you are installing a package from github and it has dependencies on Bioconductor you can use directly `BiocManager::install("user/package")` and it will use `remotes::install_github` internally and pull also packages from Bioconductor (from the right release).

BiocViews is/was used to detect that a package is for Bioconductor by remotes, devtools and sessioninfo among other packages that help the development and testing of packages. In addition Bioconductor doesn't add a `Repository: Bioconductor` on the provided DESCRIPTION field of the packages, while CRAN adds `Repository: CRAN`.

---
Tags: r, bioconductor
---


title: Mapping NCBI taxID to divergence time?

text:
# Question
Title: Mapping NCBI taxID to divergence time?
I want to map NCBI taxID, specifically a pair of Genbank sequences, to the divergence time estimate of their common ancestor. For example, using a coalescence theory approach predominantly used in phylogenetics. Do we have a database for such correlation?

# Answer
> I want to map NCBI taxID to their evolution time. Do we have a database for such correlation?

It would be best to edit your question to clarify what you mean by "evolution time" since this is not a well-defined term in molecular evolution. Even so, I think we can provide a clear answer to this question.

Estimating branch lengths between divergent taxa is an area of active research—there is no *one* method that works best in all scenarios. As far as taxonomy is concerned, resources like the NCBI Taxonomy database are under continuous refinement by the community.

Computing estimated divergence times and mapping those times to taxids would be a non-trivial task even for a single clade. Actively maintaining a database of divergence times across the entire taxonomy would be incredibly challenging, and it's doubtful that any such resource exists.
> 3 votes

# Answer
Frontiers in Ecology and Evolution did a focused and succinct review looking at 'the Tree of Life' in 2017 here. They particularly highlight the importance of integrating a fossil record within molecular trees. The upshot is simply that increasingly tighter integration between palaeontology and molecular dating is the ongoing research frontier.

The best approach in context is to look at relevant literature. The key challenge in my opinion is the underlying population dynamic because this does heavily affect temporal estimates.
> 1 votes

---
Tags: phylogenetics, phylogeny, taxonomy
---


# Question
Title: Is the algorithm of the STAR RNA Seq-Aligner similar to the Knuth-Morris-Pratt string matching algorithm?
I don't know if I'm in the right SE, if not let me know and I'll delete it.

I am reading the original publication of the STAR RNA-Sequence aligner.

Given the quoted text below I wonder if the algorithm for the seed search is actually an implementation of th Knut-Morris-Pratt algorithm?

> The central idea of the STAR seed finding phase is the sequential search for a Maximal Mappable Prefix (MMP). MMP is similar to the Maximal Exact (Unique) Match concept used by the large-scale genome alignment tools Mummer (Delcher et al., 1999, 2002; Kurtz et al.) and MAUVE (Darling et al., 2004, 2010). Given a read sequence R, read location i and a reference genome sequence G, the MMP(R,i,G) is defined as the longest substring (Ri, Ri+1, … , Ri+MML−1) that matches exactly one or more substrings of G, where MML is the maximum mappable length. We will explain this concept using a simple example of a read that contains a single splice junction and no mismatches (Fig. 1a). In the first step, the algorithm finds the MMP starting from the first base of the read. Because the read in this example comprises a splice junction, it cannot be mapped contiguously to the genome, and thus the first seed will be mapped to a donor splice site. Next, the MMP search is repeated for the unmapped portion of the read, which, in this case, will be mapped to an acceptor splice site. Note that this **sequential application of MMP search only to the unmapped portions of the read** makes the STAR algorithm extremely fast and distinguishes it from Mummer and MAUVE, which find all possible Maximal Exact Matches.

It is not cited in the references but it might be considered 'common knowledge' by now which is why the citation is omitted...

# Answer
> 5 votes
No, STAR isn't using the KMP algorithm or a modification of it. The KMP algorithm is an online exact pattern matching algorithm. It does (linear time) pre-processing on the *query* and then finds all exact occurrences of the query in the string (as you may be alluding, you can also modify KMP to report maximal exactly matching prefixes of the query). This process takes time of order N+M where N is the length of the reference and M the length of the query.

On the other hand, STAR pre-computes a suffix array on the *reference*, and then uses this to answer exact prefix match queries across *many* different queries. This is a very different approach. Building the suffix array takes time linear in the reference, and storing the suffix array requires space linear in the reference (though the constant factor is pretty big). However, it allows enumerating all occurrences of a query *Q* in the reference *R* in time $O(k + lg(\left|R\right|)+ \left|Q\right|)$ where $k$ is the number of occurences of *Q* in *R*. This holds not just for the whole pattern, but for any maximal prefix as well. Also, in practice, STAR does this much faster by building a lookup table that maps small prefixes (say, all length 12 prefixes) to the contiguous suffix array interval where they occur. In practice, this gives very fast lookup of queries. Again, a *key* difference between this kind of approach and KMP is that the suffix array pre-processes the *reference* and therefore builds an index that works for many different queries (the sequencing reads), while KMP pre-processes the query to find all of its occurrences on the reference. Another *key* difference is that, while the KMP algorithm is always linear in the length of the reference to return all matches in the worst case, the suffix array can be considerably faster.

In the original STAR publication, this is mentioned in the supplementary material in section 1.2 "Pre-indexing of suffix arrays". This idea (along with others) is also covered in quite some depth in this paper. Basically, the look up table has been a folk technique for a long time, but it's not clear there is an "original" source for it.

---
Tags: algorithms, star
---


# Question
Title: How to input data and metadata from NCBI for RNA-Seq analysis in R
So, I have downloaded some .txt and .csv files containing the raw count data (I am supposed to find out differentially expressed genes and see which genes and up- and down-regulated). I have also downloaded both the .family.soft and the series matrix files as they mentioned they contained metadata - however I'm facing several errors while trying to import them.

And all the tutorials I am trying to follow seem to have premade metadata files and there is no explanation on how this metadata is made.

Please help me in trying to understand what I am doing wrong and how to proceed. I will try to give as much information as I can to the best of my abilities, as I have barely any experience in this field.

# Answer
A metadata file isn't that large or complicated. If you want to do DE, you must have two treatments? Or two timepoints? Two or more something, and you are comparing samples of one type to another?

It's just a data table where the rownames are the samples, and then a column for, say, treatment. If you have some samples treated, some not, some at day 0, and some at day 5, then you have a column for Day and another column for Treatment. (Depending on what you want to compare to what, you might want another column of Day\_Treatment)
> 1 votes

# Answer
I don't expect the soft or matrix files will be useful for you, although it might be possible to pull metadata out of those files if absolutely necessary. I usually only see those for microarray data, so it might be worth double-checking that you do actually have RNASeq data.

I've found the best "quick start" explanation on how to carry out differential expression analysis to be the one for DESeq2. As input, this requires three things:

1. A count matrix, with samples in columns, genes in rows, and matrix entries being the raw read counts of reads mapped to genes for each sample
2. A metadata data frame, with rows arranged exactly as in the columns of the count matrix, and fields describing features of each sample
3. A statistical model / equation that describes the relationship between samples

I have found the equation to be the most complicated to understand, and unfortunately it's very situation dependent. In its simplest form, setting up a DESeq2 data structure would look like this:

```
dds <- DESeqDataSetFromMatrix(countData = count.mat,
                              colData = metadata.df,
                              design= ~ Condition)

```

This would be for a situation where all samples were sequenced in the same sequencing run, and only a single characteristic were being compared (e.g. with vs without treatment). The metadata data frame might look something like this:

```
metadata.df <- data.frame(row.names = colnames(count.mat),
                          Condition = c(rep("treated", 6),
                                        rep("untreated", 6))

```

e.g.:

```
          Condition
Sample01    treated
Sample02    treated
Sample03    treated
Sample04    treated
Sample05    treated
Sample06    treated
Sample07  untreated
Sample08  untreated
Sample09  untreated
Sample010 untreated
Sample011 untreated
Sample012 untreated

```
> 1 votes

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


# Question
Title: How to interpret the Copy Number from pseudocount depmap portal file?
*This question has also been asked on Biostars*

How do I convert the Copy Number (CN) file into useful information such as HOMDEL and AMP?

This is an already normalized file with pseudo count + 1. In this file there are values from 0 to 1 for each cancer cell line, in wide format.

I have already looked at the forum post about the file, but it's not helpful for me.

This is what the file looks like:

```
,DDX11L1 (100287102),WASH7P (653635),FAM138A (645520),OR4F5 (79501)
ACH-000621,1.0166833616923812,1.0166833616923812,1.0166833616923812,1.0166833616923812
ACH-002512,2.9659871795507433,2.9659871795507433,2.9659871795507433,2.9659871795507433
ACH-001188,0.9881267634909964,0.9881267634909964,0.9881267634909964,0.9881267634909964
ACH-000471,0.8053837915961594,0.8053837915961594,0.8053837915961594,0.8053837915961594

```

Thank you for your help

# Answer
> 1 votes
I'm not familiar with this file structure, but have you looked at the forum post about relative copy number linked to in the FAQ?

> Since we do not have matched normals, the output is a “copy ratio” or relative copy number. It is relative to the rest of the genome for that cell line. E.g. if the cell line is tetraploid we would not be able to see it from the relative copy number. These values are reported as log2(relative CN + 1) in the portal.

---
Tags: cancer, cell-line, copy-number
---


# Question
Title: How to extract and visualize sequence logos from CNN kernels?
I created a basic CNN architecture using Tensorflow to classify transcription factor binding sites. My aim is to somehow extract and visualize sequence logos from the convolutional kernels. The model looks like this:

```
model = tf.keras.models.Sequential()
model.add(tf.keras.layers.Conv2D(128, (1, 24), padding='same', input_shape=(1, 101, 4), activation='relu'))
model.add(tf.keras.layers.GlobalMaxPool2D())
# ...

```

I am using the supplement dataset from this publication

What I tried so far is that I have trained the model and created a new model out of it which has only an input layer and an output layer that contains the weights of the original convolution kernels:

```
test_model = tf.keras.models.Model(inputs=model.inputs, outputs=model.layers[0].output)
feature_maps = test_model.predict(x)

```

I have ran a prediction on a sequence and plotted out a convolution kernel:

At this point I got a bit stuck as I am not sure how can I move forward with this in the direction of visualizing the sequence logos.

# Answer
> 0 votes
What you are missing is the accuracy of the prediction, this is the next step and must be \>80%, preferably \>85% and \>90% is good in biological systems. If its lower than 80% you'll need to adjust the training. Following this there are issues with false positives and false negatives.

The method of using seems fine for mapping sequence motifs to transcription factors, which is my interpretation of your question. There must have been vectorisation of the sequence data and this is not trivial.

To answer your question, which you might be asking, you cannot produce an automated heatmap of sequence motif in relation to the transcriptor factor binding. You simply input a sequence and the associated transcription factors are revealed. Unless you vectorised your transcription factors somehow and associated this with a sequence motif.

To produce a heatmap, which again is what I suspect you are requesting, you do this via iteration. So you need a bit of coding. Thus you take all relevant sequence motifs - particularly those that have not been trained - and input them into deep learning algorithm and map which motifs associate with which binding factors and assign a numerical score to given groups of sequence motifs. Any graphics program can then be used to convert this to a heat map. Deciding 'taxonomic groups' for the sequence motifs is subject to your discretion. If you've given the transcription factors and trained this onto the sequence motifs, then classification of transcription factors I suspect would be much easier.

Heatmaps are produced in other machine learning approaches, i.e. which exclude deep learning, but your next step is percentage accuracy regardless of approach.

---
Tags: python, machine-learning, tensorflow
---


# Question
Title: Acquiring & Processing single cell data (not 10x)
I have a question related to processing of Single cell RNA-seq data originating from non 10x platform. I have worked with data that originated from 10x platform and I parse it through cell ranger pipeline and then put it through Seurat.

As an example the dataset in this GSE156644 repository are not being processed through Cell ranger pipeline and I get an error that indicates "single-end read found" another example is of this GSE115235, in the later I have bunch of files that look like bulk RNA-seq.

If you check the GSE115235 it has bunch of files that are not cellranger compatible. However, they are from single cell experiment. I find it strange that NCBI accepts such submission and the reviewers don't even pay attention to what has been submitted.

I guess, that these are SE RNA seq data that can be aligned using STAR or Bowtie and then merged in some way to get per cell per sample expression.

Along similar lines, how do I deal with 2 fastq files for a scRNA seq run? example https://www.ebi.ac.uk/ena/browser/view/PRJNA602526

This link has scRNA seq files that look paired end (R1 & R2 but NO I1). I am bit confused as cell ranger count option requires the folder to have 3 files (I1,R1,R2). How can I run cell ranger with R1 & R2? Or do I need to use another aligner?

I will send an email to sra and also to the authors, but any other help to navigate through these datasets and process them is deeply appreciated.

Thank you.

# Answer
> 2 votes
*Modified from @devon-ryan's comments:*

The couple of samples I looked at have the sequence from read 1 embedded in their header. I1 isn't needed, they're already demultiplexed. You could feed those into STAR or similar tools, but processing them will be more of a pain since most tools expect the original unprocessed version of what they uploaded.

*Modified from atpoint's comments:*

Single-cell data is not really useful as "bulk" as the notable amount of PCR during the library prep creates quite a bias that without the UMI information cannot be corrected. I'd email the authors and ask for the R1 files. Single-cell uploads at NCBI tend to be a mess, that's unfortunately how it is. Without R1 data are almost useless.

Those file format issues are strange, but it is what it is. Reviewers do not pay attention to method sections, choice of tools and data submission most of the time in my experience, especially not when the paper is biology-focused rather than a computational-focused one. Even more, you often do not find the suitable metadata for a proper analysis such as information about batches etc. The only choice you have is emailing the authors as noted above.

If by "two files" you mean that the same library (=sample) was sequenced over multiple lanes (very common to do that since required depth is often larger than what a single lane can provide) then merge the fastq files, so all R1s into one and all R2 into one file.

---
Tags: scrnaseq
---


# Question
Title: How to improve a genome assembly using Dovetail and PacBio assembly?
I have more of a conceptual question. I have two genome assemblies from the same plant, one from Dovetail technology (~998 Gb) and another is PacBio HiFi assembly (~1.1 Gb). The Dovetail assembly is more contiguous but has lower base quality whereas the PacBio has higher base quality but is more fragmented. It is a diploid organism with relatively high heterozygosity.

I tried to use `RagTag` `patch` to improve the assembly. But no improvement was obtained.

Is there a way to use both assemblies and produce a hybrid assembly with high contiguity and base quality?

# Answer
> 1 votes
*Answer from @maximilian-press, converted from comment:*

I am confused when you talk about base quality, because those technologies are largely about ordering and orienting existing contigs, not assembling contigs in the first place. If you have Hi-C reads, you can read them into juicebox and actually visualize the junctions to see if the assembly makes sense: aidenlab.org/juicebox. It would also be helpful to see the actual numbers, not just "more" or "lower".

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


# Question
Title: Plink error: --check-sex/--impute-sex requires at least one polymorphic X chromosome locus
I am running an analysis in plink.

This is my first attempt at a QC and I keep encountering the same message when I try to go through with my sex check. I've done the SNP missingness step, it's this step that keeps stumping me. I'm assuming it has something to do with the chromosome naming but I'm a bit lost as to how I should tackle that problem. For reference, I attempted running the --check-sex function on the files (which are in bim, bed, fam) however, I'm met with this error.

```
431911 variants loaded from .bim file.
13635 people (7170 males, 6465 females) loaded from .fam.
13635 phenotype values loaded from .fam.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 13635 founders and 0 nonfounders present.
Calculating allele frequencies... done.
Total genotyping rate is 0.997824.
431911 variants and 13635 people pass filters and QC.
Among remaining phenotypes, 4457 are cases and 9178 are controls.
**Error: --check-sex/--impute-sex requires at least one polymorphic X chromosome
locus.**

```

Does anybody know what I should do to fix this?

This is the command I used: plink --bfile ppmissingness2 --check-sex

# Answer
Without X chromosome variation, sex cannot be reliably determined from input genotypes.

The error is telling you that there are no polymorphic X chromosome loci in the provided data. Either your input files don't include any X chromosome information (this could happen due to the increased statistical challenges involved in X chromosome data), or all the X chromosome "variants" are invariant / homozygous.
> 1 votes

---
Tags: gwas, plink, bed, qc
---


# Question
Title: What does it mean by "PC" in this article?
What does it mean by "PC" in this article?

# Answer
*Answer from @user438383, @matteo-ferla and @m, converted from comments \[anyone, feel free to edit this answer to improve it\]:*

According to the intro it means 'pseudocenter'. All measurements are relative to an assigned point because a protein doesn't have a 'centre' they assigned it as a pseudocentre.

Note that PC normally means principal component, but that does not seem to be the case here.
> 3 votes

---
Tags: protein-structure, algorithms
---