Spliced vs. unspliced ratios for transcripts in RNA-seq data

I don't know about tools, but i've used the following python code to calculate the ratio of reads that overlap the 5' or 3' ends of introns or that are spliced. We sum these across all introns in a gene set (we actaully use this for iCLIP analysis to see if RNA binding proteins bind pre-mRNA or spliced RNA).

import pysam
from collections import Counter
from CGAT import GTF, IOTools

def calculateSplicingIndex(bamfile, gtffile, outfile):

    bamfile = pysam.AlignmentFile(bamfile)

    counts = Counter()

    for transcript in GTF.transcript_iterator(
            GTF.iterator(IOTools.openFile(gtffile))):

        introns = GTF.toIntronIntervals(transcript)

        for intron in introns:
            reads = bamfile.fetch(
                reference=transcript[0].contig,
                start=intron[0], end=intron[1])

            for read in reads:
                if 'N' in read.cigarstring:
                    blocks = read.get_blocks()
                    starts, ends = zip(*blocks)
                    if intron[0] in ends and intron[1] in starts:
                        counts["Exon_Exon"] += 1
                    else:
                        counts["spliced_uncounted"] += 1
                elif (read.reference_start <= intron[0] - 3
                      and read.reference_end >= intron[0] + 3):
                    if transcript[0].strand == "+":
                        counts["Exon_Intron"] += 1
                    else:
                        counts["Intron_Exon"] += 1
                elif (read.reference_start <= intron[1] - 3
                      and read.reference_end >= intron[1] + 3):
                    if transcript[0].strand == "+":
                        counts["Intron_Exon"] += 1
                    else:
                        counts["Exon_Intron"] += 1
                else:
                    counts["unspliced_uncounted"] += 1

    header = ["Exon_Exon",
              "Exon_Intron",
              "Intron_Exon",
              "spliced_uncounted",
              "unspliced_uncounted"]

    with IOTools.openFile(outfile, "w") as outf:

        outf.write("\t".join(header)+"\n")
        outf.write("\t".join(map(str, [counts[col] for col in header]))
                   + "\n")

Unfortunately this uses a bunch of libraries you may or may not have, including CGAT (for the GTF parser and IOTools package) and pysam.

Once you've got these statistics you can calculate the "splicing index" as
the log2 ratio of 2 times the number of spliced reads divided by the number of reads overlapping the 3' and 5' ends of introns.


Python: How to write duplicate sequences removed from fasta file to new file

Please do not reinvent the wheel - seqkit rmdup does what you want.
From the website:

rmdup
Usage
remove duplicated sequences by id/name/sequence

Usage:
  seqkit rmdup [flags]

Flags:
  -n, --by-name                by full name instead of just id
  -s, --by-seq                 by seq
  -D, --dup-num-file string    file to save number and list of duplicated seqs
  -d, --dup-seqs-file string   file to save duplicated seqs
  -h, --help                   help for rmdup
  -i, --ignore-case            ignore case

How does a -D file look like, you ask:
$ zcat hairpin.fa.gz \
    | seqkit rmdup -s -i -o clean.fa.gz -d duplicated.fa.gz -D duplicated.detail.txt

$ cat duplicated.detail.txt   # here is not the entire list
3   hsa-mir-424, mml-mir-424, ppy-mir-424
3   hsa-mir-342, mml-mir-342, ppy-mir-342
2   ngi-mir-932, nlo-mir-932
2   ssc-mir-9784-1, ssc-mir-9784-2


Ironically, I found that this problem was caused by the fact that the tables are too large!
They really should have designed the correct error messages!


How do you generate read-length vs read-quality plot for long-read sequencing data (e.g., MinION)?

I also wrote a package to create various plots from Oxford Nanopore sequencing data and alignments: NanoPlot. It can be installed through pip (see also the README on Github). In addition to multiple plots also a limited NanoStats output is created (see also NanoStat). Data can be presented using:

A fastq file (optionally compressed)
A bam file
The sequencing_summary.txt file generated by albacore

Using optional flags you can:

Log transform the read lengths
Use aligned reads rather than sequenced reads
Downsample the reads
Set a maximum read length

I've added an example below, plotting log transformed read length vs average read quality (using a kernel density estimate). More examples can be found in the gallery on my blog.

I welcome all feedback and suggestions!


Split fasta file based on groups in header information and output as separate files

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

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


Have a read of the Estimation of model accuracy in CASP13 paper to get the latest rankings.

Generally you would run the model quality software on the models you have generated and choose the best. Also have a look at the structures in PyMol, you learn a lot that way. See, for instance, how much variety there is, and where the models vary. Is there a conserved core in many of the models? Or is there just a blob of "spaghetti" that doesn't look like an ordered protein?

If you have many models (e.g. from a Rosetta run) you might want to look into clustering the models first to reduce the number you have to work with.


How to initiate RNA-Seq analysis of TCGA files?

The GDC API user guide (https://docs.gdc.cancer.gov/API/PDF/API_UG.pdf) provides some guidance on the structure of GDC data. So for getting RNA-Seq data from GDC, you would provide this filter to get the 551 text files with workflow_type set as "HTSeq - FPKM-UQ" and "TCGA-LUSC" for the project. There are other workflow types as well:

cases.project.program.name in ["TCGA"]
and cases.project.project_id in ["TCGA-LUSC"]
and files.access in ["open"]
and files.analysis.workflow_type in ["HTSeq - FPKM-UQ"]
and files.data_format in ["txt"]
and files.data_type in ["Gene Expression Quantification"]
and files.experimental_strategy in ["RNA-Seq"]

(You don’t need all those filters to get the 551 files; I included them to show what is available.) The txt files are returned as gzipped tar files that hold text files with rows of the form:

ENSG00000158486.12  5341.21601153

The above filters provide all the FPKM-UQ files for TCGA-LUSC. If you wanted to just get the data for normal samples as controls, that would require filtering on sample type as well, adding in the additional filter:

and cases.samples.sample_type in [“solid tissue normal”]

That will return 49 files. If you instead wished to get only the tumor files, the additional filter would instead be:

and cases.samples.sample_type in [“primary tumor”]

and that returns 502 files. Note that there are several different sample types for tumor and normal tissues; those just happen to be the ones used here.
(Note that if you are using the GDC Repository portal, adding that additional filter requires you to select “Add a Case/Biospecimen Filter” on the “Cases” filter tab.)

Since GDC is part of the NCI Cancer Research Data Commons (CRDC), instead of downloading the data, you can actually work with it in the cloud using the various CRDC Cloud Resources such as the ISB-CGC: https://www.isb-cgc.org/. In fact, all the RNA-Seq data for TCGA is in an ISB-CGC Google BigQuery table: (https://console.cloud.google.com/bigquery?p=isb-cgc&d=TCGA_hg38_data_v0&t=RNAseq_Gene_Expression&page=table)
And the above count data can be pulled using the query:

SELECT HTSeq__FPKM_UQ
FROM `isb-cgc.TCGA_hg38_data_v0.RNAseq_Gene_Expression`
WHERE Ensembl_gene_id_v = "ENSG00000158486.12" AND project_short_name = "TCGA-LUSC"

If you want to discriminate between normal and tumor expression levels, the ISB-CGC BQ tables make it pretty easy to tag the expression values with the type of tissue involved. The sample type can be assigned by joining the previous query with the biospecimen table:

WITH bio as (
SELECT sample_barcode, sample_type_name
FROM `isb-cgc.TCGA_bioclin_v0.Biospecimen`
),
expr as (SELECT sample_barcode, HTSeq__FPKM_UQ
FROM `isb-cgc.TCGA_hg38_data_v0.RNAseq_Gene_Expression`
WHERE Ensembl_gene_id_v = "ENSG00000158486.12" AND project_short_name = "TCGA-LUSC")
SELECT bio.sample_type_name, expr.HTSeq__FPKM_UQ
FROM expr
JOIN bio
ON bio.sample_barcode = expr.sample_barcode

The GDC has listed some of the available sample types on this page: https://gdc.cancer.gov/resources-tcga-users/tcga-code-tables/sample-type-codes
Also, BigQuery has a user interface that makes exploring data sets and writing queries easy, but you need to log in with a Google ID, so the ISB-CGC has a BigQuery Table Search (https://isb-cgc.appspot.com/bq_meta_search/) that allows you to explore the available tables without logging in. If you want to run queries, Google offers a free tier that you can use to try things out and run small jobs (https://cloud.google.com/free/) and ISB-CGC offers free credits too. ISB-CGC also has examples of working in R with BigQuery within the Community Notebook Repository (https://isb-cancer-genomics-cloud.readthedocs.io/en/latest/sections/HowTos.html). A few useful notebooks would be 'How do I create cohorts of patients?', 'How do I plot a heatmap using data in BigQuery?', and 'How do I get started fast?'.


How to tell if library generation for RNA-seq experiments are stranded vs not stranded

If you already have a reference transcriptome assembly, then you may run the Salmon tool for quantification using your RNA-seq samples. It has an automated library type inference function (https://salmon.readthedocs.io/en/latest/salmon.html#what-s-this-libtype). It might help you to determine whether the libs are stranded or not.


How to get 2d/3d structural coordinates from a PDB file?

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.


Biopython: resseq doesn't match pdb file

Let's check it. I copied your file fragment into a.pdb.

>>> from Bio import PDB
>>> for res in PDB.PDBParser().get_structure('pdb', 'a.pdb')[0].get_residues():
...   print(res.get_full_id())
...
('pdb', 0, 'A', (' ', 22, ' '))
('pdb', 0, 'A', (' ', 23, ' '))

So BioPython (ver. 1.66) gives 22 and 23 as expected.


In matlab this can be done using the SimBiology toolbox.
In mathematica this can be done using the SystemModeler.
In python there are multiple packages, e.g. (in no particular order) PySB, Tellurium, PySD


The easiest way is to check annotations in mmCIF or PDB files.
Looking at my page with statistics of mmCIF tags, last updated in March:
_struct_conn_type.id        59137×metalc 42706×covale 36195×disulf 10725×hydrog
_struct_conn.conn_type_id   5350402×hydrog 1733277×metalc 559092×covale 226472×disulf

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


alternatives to PHASTER phage annotation server

The answer is to use phastaf, a recently developed software tool developed by Torsten Seemann (the same developer as prokka), which comes with a packaged phaster.ca db.  Get it here


You can use the "manual network construction" approach. Calculate TOM, then use exportNetworkToCytoscape or exportNetworkToVisANT to export the network (or simply save the TOM matrix as a table and load it in Python).


Bulk RNA-Seq Read Length Normalization across different samples

The simplest thing to do is to trim the 150 bp fasts so that they are 100 long.  I don't think there is an easy way to correct for the fact that the 150 bp long reads will have a higher unambiguous rate of alignment and gene assignment than the 100 bp long reads.
If you have a mix of all the experimental conditions across all the lengths, you can include length as an element in the design, so that DESeq2 may be able to model the difference between the two batches.


Differentiating molecules based on peptide sequence? How to annotate?

Knowing almost nothing about the biology, I think that the most straightforward way to do this is:

train profile HMMs on your different classes of protein sequences (in human maybe, but ideally including MHCs from a large number of organisms)
translate your transcripts into proteins and run them against your trained profile HMMs
confirm that the profiles adequately discriminate the two classes and adjust if needed.

For an example of how to do this (on MHC), you can see here.
These approaches are then easily extensible to further organisms, as you have the profile HMMs.
It's possible that a simpler approach would consist of just finding hits against human MHC and then just checking for the existence of sequence patterns corresponding to those structural motifs, but I doubt that would work well in any significantly divergent organism.


GMAJ error: Exception in thread "main" java.lang.ClassCastException: java.lang.StringBuffer cannot be cast to java.lang.String

It's a programming error that should look like this:

StringBuffer sb = (...)
(...)
String s = (String)sb;

while it should be

StringBuffer sb = (...)
(...)
String s = sb.toString();

( see also https://www.biostars.org/p/177992/ ) . You cannot do anything about it but asking the authors of the software.


What is "Possible polymerase run-on fragment"?

On most coding genes (with the exception of replication dependent histone genes), Transcripts (as opposed to transcription) are terminated by cleavage polyadenylation sites, where the growing transcript is cleaved by a protein complex and the poly-A tail. However, PolII continues to transcribe until it is terminated by torpedo termination (see here). Normally this extra RNA after the polyA site is degraded by Xrn2, but if you capture any particular snapshot of the RNA in the cell at any one time there will be a small amount that corresponds to this run on transcription.


Duplicated genes and genome fragmentation

The tBlastn looks good and looks like some sort of parology, although the hits you are getting need to be checked to make sure they are not psuedo-genes. The shot-gun genome however appears poorly annotated because it had misannotated 5-1 and 5-2. What I don't understand is why you haven't performed a blast back over the original highly annotated genome? You must have done this?

There are a few avenues to go down here, what I'd like to know which group of protists are you working with? I don't need to know the species, but I do need to know the family please.


Using `install.packages` with conda-managed R

Installing packages outside of an R session is arguably the point of using conda.  Package management isn't something that should be done within scripts, only the loading of packages. Doing so ensures that your environment setup is somewhat robust, which should be set before you perform actual calculations.

Using conda install ... is the best way to ensure environment stability here, and, thankfully, lots of people at Bioconductor have worked very hard to make sure bioconda packages fit with the conda framework.

If you need to find packages, you can search online at https://anaconda.org, or from the command line via conda search <pkg>.


Pipeline to analyse small RNAseq and tRNAs?

I have heard that the Small RNA Workbench from UEA can be an easy solution for a full small RNA pipeline.  It is java based, has a GUI, and is supposed to be easy to install and run.  I am not an expert in small RNA analysis, but I think this pipeline will certainly given you an idea about the steps needed for a small RNA analysis, it has been put together by experts.


How can I run the latest version of samtools on Slurm?

I would start with:
module avail

Then look to see (a) if there is a samtools module avilable and (b) exactly what it's called.
If it's called samtools/1.11 then:
module load samtools/1.11

Should work, if it doesn't speak to your sys admin.


Plotting metadata from GenomicRanges object

You can access the data.frame just using $. All of them can be accessed using values() . Something like this:
library(GenomicRanges)
gr = GRanges(seqnames=1,IRanges(start=1:10,width=1))
values(gr) = data.frame(eigen=rnorm(10),ratio=runif(10))

plot(start(gr),gr$eigen)

or
plot(start(gr),values(gr)$eigen)


How to generate Genbank format out of paired end fastq data?

The steps needed to assemble a genome in Galaxy are layed out in the Assembly training material from the Galaxy Training network. The software for this should be available on the main galaxy servers (e.g., usegalaxy.eu). Since you have a prokaryote, I suggest you follow the "Unicycler Assembly" workflow, which is also available on usegalaxy.eu.


Separating peaks of chip-seq with specific length

You can do it by awk.
awk -v OFS="\t" '{ if ( ($3-$2) >= 500 && ($3-$2) <= 1000) { print } }' in.bed > out.bed

Assuming there is no header.


GATK documentation for required depth to reliably call heterozygous mutation in diploid organism?

UPDATE: As an update, Sarah Walker (co-author on the poster) responded to my question on the GATK forum. She clarified with the following statement:

  We believe the sites around 30X (and above 150X) are being filtered
  due to low mapping quality (since it is whole genome so there are many
  areas that are hard to map), which explains the low sensitivities in
  these areas. As the error bars show, there are very few variants in
  these low mapping quality areas.

She then posted the following plot that falls much more inline with what we would all expect. This is a plot showing sensitivity for NA12878 compared to NIST truth.

shlee, from the GATK forums, pointed me to a poster from AGBT 2018 where the authors compared the HiSeqX and NovaSeq instruments. Figure 3 shows sensitivity for SNPs and short INDELs across multiple allelic fractions (AF) from spike-ins. I was surprised to find that SNPs where AF = 0.4 only had ~50% sensitivity at ~25x coverage, and only ~75% sensitivity at ~50x coverage.

The poster is targeted at somatic variation, but data from figure 3 should behave like a standard diploid germline variant because you often see a wide range of allelic fractions from germline variants. They used spike-ins for figure 3. It's not clear to me whether they used HaplotypeCaller in figure 3, but I don't see how they could have used Mutect2, like they did in figures 1 & 2.

Based on this information, it appears the standard 30x genome-wide coverage is generally insufficient to confidently call many heterozygous mutations, and that ~90-100x is ideal.

Am I missing something? Are we really missing that many heterozygous mutations in standard WGS/WES?

I'll leave this question unanswered for now to allow others to chime in on interpretation, etc. I'd appreciate additional thoughts.


FYI coronaviruses have proof reading activity, unlike HIV, its a huge genome, unlike HIV. Thus the point mutation rate is unlikely to be as rapid as HIV.

The phylogenies in this paper are not really about the evolution or epidemiology of the virus. The trees used for viral taxonomy and form part of the International Committee for the Taxonomy of Viruses (ICTV) assessment in officially naming the virus. The ICTV is the authority in naming a virus. The ICTV essentially called it SARS2 as opposed to the WHO who used the name COVID-19. Thats the point of this paper and will be massively cited.

The broader evolutionary epidemiology of "SARS2" is for further studies. The Lancet paper on 22nd February gave some indication of this (and I suspect you know about that paper).


Interpreting bayesian phylogenetic model using Tracer?

When I have a moment I will fill in the blanks

You must have an RNA virus (nothing else will mutate that fast). The basal group and root is the USA. It estimates the time to that basal event is 32.3 (likely years).

"frequencies" frequency of ACGT where the missing one will be 0.25.

Its a simple model of constant population size.

ESS is the most important parameter to ensure convergence of the Bayesian simulation .

Likelihood is the second most important parameter. I'll edit again shortly.

I'll fill in one more blank and summarise ...

Kappa ...
This will be the transition transversion ratio in HKY. I suspect it also will be transition transversion ratio in Kimuras 3-parameter model, this is because ... the transversions are further split into two mutation groups, ie. purine to pyrimidine has 2 further combinations and you end up with 2 ratios.  that makes sense, transition/transversion (kappa 1), transition/transversion 1 (kappa 2), transition/transversion2 (kappa 3) ....

CP1.kappa 18.716
CP2.kappa 17.16
CP3.kappa 22.796

CP2 and CP3 are either side of CP1 .. which is what you would expect on my interpretation.

These are very high values by the way, essentially there are very few transversions in this data, so ... I dunno. You need to check the data, if you see lots (and lots) of A-G and T-C mutations, then thats what I believe it is saying.

Briefly posterior and prior are the fundamental method of a Bayesian calculations, the output is probability expressed as a log, so the more negative the higher the probability (the prior looks like a random tree, or a reasonable tree with really rubbish parameters, likely the former).

Summary
Essentially I am going to call it a day here. The real issue is that the output summarise the entire modern theory in phylogenetics, so its a very long response required to address all of the parameters and I honestly don't think it helps your analysis. What I am clear about is:

* You've omitted a gamma distribution (I think Tracer would summarise this) *

So my advice is:

ESS MUST be homogeneous between runs - big trouble if not
gamma distribution MUST be run

Beast can be tricky getting gamma to converge via ESS, just a word of warning.

I can't see the full robustness, because its absent... However, from the output, I think its pretty clear, that you (well me) can more or less see the model you ran.


Tool to remove a PCR contamination in NGS data

This bash script will get you all of the reads that do not start with your primer sequence. It also handles degeneracy. It requires the tool bioawk.

It works by looking at the first bases in the read and only prints out the headers that do not have a match with the degenerate primer. Then, comm makes sure that the read was retained in both the R1 and R2 filtering. This allows reads to be filtered whether or not the primers start at the beginning of the R1 or R2 read.

seqtk subseq is where you recover the reads that don't start with the primer.

Be warned though, any off-by-one or sequencing errors will allow some of the sequences with primers at the beginning through. Also, this script requires that the headers for the R1 and R2 reads are the same and aren't appended with anything like "\1", "\2", et cetera.

#!/bin/bash
get_primer_free_starts () {
    bioawk -cfastx 'BEGIN{seqs[0] = "CATGGGAGGGCGATCGCAA"; seqs[1] = "CCTGGGAGGGCGATCGCAA"; \
        seqs[2] = "CGTGGGAGGGCGATCGCAA"; seqs[3] = "CTTGGGAGGGCGATCGCAA"} \
        {thisseq = substr($seq, 1, length(seqs[0])); printme=1\
    for (i=0; i<=3; i++) { \
        if (thisseq != seqs[i]){ \
            printme=0 \
        }
    }; if (printme==1){ \
        print($seq)
    }}' "${INPUT} | sort > "${OUTPUT}
}

INPUT="R1.fastq.gz"
OUTPUT="R1_OK.txt"
get_primer_free_starts

INPUT="R2.fastq.gz"
OUTPUT="R2_OK.txt"
get_primer_free_starts

comm -12 R1_OK.txt R2_OK.txt | sort | uniq > passing_reads.txt

seqtk subseq R1.fastq.gz passing_reads.txt > R1.filtered.fq
seqtk subseq R2.fastq.gz passing_reads.txt > R2.filtered.fq


As was suggested by @benn, you will need to add the TdTomato sequence to the end of your fasta formatted genome file. If your TdTomato sequence is in fasta format already, cat will do the trick. If not, just go the end of the file, start a new line and add a header line starting with ">" then immediately the name of your gene, no whitespcaes in between, something like ">TdTomato". Then add your sequence in a new line.

The format of the GTF file is explained in the image you posted, the trick is to use the exact same name (in column 1) that is used in the genome file you will create at the previous step (no ">" though). Again, this line, in which the different fields (columns) will be separated by tabs, will be added to the end of your GTF file.


How to use SBOL (Synthetic Biology Open Language)?

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.


Variant calling for a subset of genes using whole gneome sequencing data

Your approach is not recommended. Always align against the entire genome. If you align against a subset you might accumulate false mappings. The aligner will always try to find the best match for every read. If the true origin of the read is not in the reference then the aligner will still try to find an acceptable mapping position. Therefore, align against the full reference and then extract the regions you are interested in, e.g. with the region option of samtools view:

samtools view [options] <in.bam>|<in.sam>|<in.cram> [region ...]

Check its manual for details.

While there is basically nothign wrong with VarScan it is quite old and not maintained anymore. I suggest you give bcftools a try (check its manual). Other maintained alternatives can be (among others) freebayes, strelka2, VarDict or the GATK from the Broad institute.


How can I downsample a BAM file while keeping both reads in pairs?

samtools has a subsampling option:

  -s FLOAT:
  Integer part is used to seed the random number generator [0]. Part after the decimal point sets the fraction of templates/pairs to subsample [no subsampling]

samtools view -bs 42.1 in.bam > subsampled.bam

will subsample 10 percent mapped reads with 42 as the seed for the random number generator.


There is now a post on Biostars:

  We are starting to define the release process for the Human Cell
  Atlas.

  The first couple of HCA releases will be primarily single-cell RNA-Seq
  and we are interested to know what you would expect our releases to
  contain

  Please fill out our survey if you would like to help us define our
  releases

  https://www.surveymonkey.com/r/S6Y753P

And I suppose that answers the original question: the data has not been released yet and probably won't be for a while considering they are still discussing how to make it available.

Update (from the same post):

  The first production release of the DCP is happening in June.

  The first data release of HCA is targetting Q4 2018.


Can an alternate allele be more common than a reference allele?

Yes, depending on how the reference allele was determined. If, for example, you want to refer to hg38 as the reference genome, then it's important to use the hg38 allele as reference.

Without stating the reference allele convention? Not that I'm aware of. However many people refer to minor/major alleles rather than non-reference alleles, in which case it is assumed that the major allele is the most common allele.

As long as you state what your reference genome and/or population is, then it's fine to choose your own convention about reference alleles. But please be consistent and explicit about what you've done. One of my biggest struggles in GWAS has been working out the correct strand orientation for stated variants (important when merging datasets), and it doesn't help for complementary variants (i.e. A/T, C/G) when the ref/alt alleles are arbitrarily chosen.


ggplot2 geom_bar() , do not merge y-values shared same x-ids

As you seem to have surmised, the graph is not erroneous, but correct. You have a total of 10 for a and that's what it's showing. If you want it to ignore duplicate X-axis values then you need to tell ggplot which datapoints it should actually be using. barplot() is assuming that your duplicate a values are intentional, ggplot assumes they should be grouped. These are both reasonable assumptions and it's your responsibility to determine beforehand which is correct.


Differential expression analysis when nested effects

llrs is right.  "Matrix Model not full rank" means that two of your conditions are confounded.  You can't get your normals from one place and your experimentals from another

And no, there generally is no clever way out of this.  For all you know, any differences you see between your normals and your cancers are 100% due to being sequenced somewhere else, and have nothing to do with real biology.

If you want to run this through DESeq, drop batch effect.  But tell whoever you are reporting the results to that batch is hopelessly confounded with cancer state, and you can't tell which differences are due to what.


Wow, I would just go for accuarcy as the singular most important measure, but be careful of sampling bias.

If the classification training had an unbiased sample* and you have an accuracy of >90%, in general ANN that wouldn't be great, but in bioinformatcs thats really good. If your training comprises sampling bias that gets tricky and you either need to accomodate that in your accuarcy ROC etc ... or else retrain to eliminate the sample bias.

You can account for all of them of course by wrapping everything into a single analysis.

*, by this I mean that each of the 3 classifications comprises 33% in each training set. If this is badly skewed, this requires consideration and gets complicated. If it is precisely 1/3 each then 'accurary' is the key stat.


In case it is still helpful as the post is rather old, the code below would generate a heatmap with annotations thanks to the ComplexHeatmap package.

But before that I would like to stress out the importance of example data, I bet a lot of experts could not help you just because of they did not have single cell data at their disposal.

I am using some scRNA-seq data stored within a Seurat object. The expression data come from the @scale.data slot and the annotations are from the @meta.data slot. I am using a few general cell type markers as genes and 3 columns from the metadata slot (cell type, cell subtype, sample region). When using the ComplexHeatmap package, the trick is to have the same order in the annotation and actual data tables and this is how you should control what (in your case single cells) is plotted in what order.

library(Seurat)
library(dplyr)
library(ComplexHeatmap)
library(circlize) # for the colorRamp2() function

seurat_object <- readRDS(...)

# Annotations

meta_data <- seurat_object@meta.data

# order of annotations/colors are defined here
ordered_meta_data <- meta_data[order(meta_data$cell_type), ]

# OPTIONAL: YOU CAN PICK COLORS FOR EACH LEVEL OF ANNOTATION
# HERE I PROVIDE COLORS FOR ONLY TWO OF THE THREE LEVELS, COLORS FOR THE
# REMAINING LEVEL IS TAKEN CARE OF BY THE PACKAGE
annotation_colors <- list("cell_type"= c("Epithelial" = "red",
                                         "Immune" = "orange",
                                         "Other" = "gray",
                                         "Stroma" = "green"),
                          "sample_region" = c("B" = "blue",
                                              "C" = "yellow",
                                              "N" = "maroon",
                                              "T" = "darkblue"))
ha = HeatmapAnnotation(df = ordered_meta_data,
                       show_annotation_name = TRUE,
                       col = annotation_colors)

# Expression data

genes_to_use <- c("EPCAM",
                  "CD3D",
                  "NCR1",
                  "CD79A",
                  "LYZ",
                  "COL1A1",
                  "CLDN5",
                  "S100B",
                  "KIT")

seurat_object <- ScaleData(seurat_object,
                           genes.use = genes_to_use)

my_data <- seurat_object@scale.data

# COLUMN ORDER OF THE EXPRESSION DATA SHOULD MATCH THE ROW ORDER OF THE
# ANNOTATION TABLE
my_data <- my_data[, rownames(ordered_meta_data)]

# Heatmap

col_fun = colorRamp2(c(-2, 0, 2), c("blue", "white", "red"))
png(filename = "example_heatmap.png",
    width = 1000,
    height = 1000)
Heatmap(
  my_data,
  col = col_fun,
  cluster_rows = FALSE,
  cluster_columns = FALSE,
  column_order = NULL,
  show_row_dend = FALSE,
  show_column_dend = FALSE,
  show_row_names = TRUE,
  show_column_names = FALSE,
  use_raster = TRUE,
  raster_device = c("png"),
  bottom_annotation = NULL,
  top_annotation = ha
)
dev.off()


Holistic enzyme activity determination with computation

Your title says holistic. This is a tad problematic as there's layers upon layers. Say, post-translation regulation, inhibiting metabolites, interacting protein etc. That is why when talking of an enzyme inhibitor, at the biochemical level one speaks in terms of k_i (inhibition constant), while at the cellular level ("holistic") one talks of IC50.
This could be a biocatalysis/synthetic biology problem (make an enzyme better) or translational medicine problem (why is this mutation pathogenic).
In terms of in silico enzyme design, one can have all sorts of issues that need to be addressed:

catalysis
side reactions
stability
aggregation
solubility

Also there's the expression problems, which may be due to subpar promoter, insufficient cofactors and poor codon usage.
For clinical variant assessment, a variant may:

destabilise leaving to degradation
remove a post-translational modification
disrupt the active site, affecting catalysis
lose a binding partner
result in incorrect localisation
result in solubility issues

So force-field scores are not the sole solution. Furthermore, a forcefield method gives one a ∆∆G (difference in relative Gibbs free energy), not a K_M (Michaelis constant) and a k_cat (catalytic rate), although it is somewhat possible to calculate these, or a relative variant of these.

a variant with just 1 or 2 variants

A variant protein can have one or two mutations (e.g. G10R) in its sequence. A mutant and a variant are near synonyms, but variant is generally intended to be a mutant in a population.

Does these application be used for such purpose?

For enzyme design MD simulations are not generally the first step. Rosetta pmut_scan for example is a nice starting step as the mutational landscape of a protein is N^19, while a MD simulation takes weeks and 6+ core per single variant.
This calculates only an increase in stability (relative difference in Gibbs folding energy), not reactivity. For that you can use a transition state, whose results are inaccurate, but majorly affected by the set-up, i.e. junk-in junk-out plays a big part. There's a saying 'a week in the lab saves you an hour at the library': this is backwards for sarcasm as it happens too often, namely that is background research is always a lifesaver. An QM-MD simulation on the wild type reaction definitely can help in elucidating things if it is not clear.
Online tools
Generally, however, there are resources out there to do many stability calculations already.
SynBio If you are interested in picking a residue or two to make an enzyme more stable, without resorting to a thermophilic enzyme, PROSS is a very nice tool.
Clinical (Caveat: If your variant is a variant that affects you, a friend or a relative, you/they should speak to a genetic councillor and not do research on it yourself) There are several tools that can describe a variant at the protein level. Do remember that ∆∆G is a potential so more negative is more stable —except in SDM where the score is inverted.

Missense3D
Venus
Miscast
SDM


Manta SV caller: I can't find a file in the path specified by the user guide

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.


Combat for multiplatform batch correction

Usually with microarrays you want to make a case/control comparison, so I am going to assume that.

Data from different array platforms is generally difficult to compare:  each platform is measuring potentially a different part of the expression of the gene (different exons or different regions), each platform is likely to need different normalisation approaches, and have different dynamic range.

If you have case/control data from several different platforms, then the most robust thing to do is to analyse each platform separately, and produce a list of up and down-regulated genes according to each platform.

Then, it would be reasonably safe to combine the lists of genes in a 'meta-analysis' approach.  I would expect that in a lot of cases, 'up' and 'down' regulation would be correlated across platforms, but I would not necessarily expect very strong correlation in the magnitude of changes between platforms for individual genes.


Is outlier detection the right choice for determining if a value lies singificantly away from the mean in a population in gene frequency?

I agree you have the right idea because that data doesn't look at all linear, so regression analysis looks really complicated. Your approach a straight comparison keeps it simple and is within the biological question you are asking,

However, I would use a non-parametric such as Wilcoxon instead of a T-test. I don't fully know the null distribution you are comparing it against. Wilcoxan doesn't require an assumption of normality. T-tests invariably require a transformation to normality and that gets complicated, but a log or sqrt transformation might be okay. If Wilcoxan produces a no significant difference you report it and move onto a T-test.

It is difficult to use outlier analysis with 9(?) observations, thats used if you have a way-out result and can afford to exclude it (>3 standard deviations). Outlier analysis is to remove a single result that is messing up your analysis and doesn't affect the bulk of the data. You alot of observations to be considering outlier analysis.


Among other things, the canon-gff3 program in the AEGeAn Toolkit infers and prints missing features like UTRs.

$ canon-gff3 < genes.gff3 > genes-with-utrs.gff3


I would ignore peak calling for this and instead compute enrichment of ChIP/input for the genome (e.g., with deepTools or presumably homer) and then plot it for the genes of interest individually (e.g., using IGV or pyGenomeTracks) or as a group (e.g., with computeMatrix). If the peaks are obvious and you trust your peak calling then sure you can just use that instead, but given that that's prone to error it's usually better to simply use normalized signal.

Regarding actually looking at occupancy, I don't think that's useful. Occupancy isn't binary. If you're instead wanting to associate peaks with genes then using the nearest gene is common (it's not great, but it's not terrible either).


SAM format: Does the BAM "Integer or numeric array" field no longer exist? Why?

An example would be ZC:B:i,271,271,1,0, which is from TMAP. The i part denotes the values held in the array, which are comma-separated.
That code is only used in samtools fastq, where that tag type isn't output into fastq files.


Tools for simulating Oxford Nanopore reads

Simulators designed specifically for Oxford Nanopore:

NanoSim
NanoSim-H
SiLiCO
ReadSim
DeepSimulator

General long read simulators:

Loresim
Loresim 2
FASTQsim
LongISLND

For an exhaustive list of existing read simulators, see page 15 of my thesis, Novel computational techniques for mapping and
classifying Next-Generation Sequencing data.


Obtaining uniquely mapped reads from BWA mem alignment

To exclude all possible multi-mapped reads from a BWA-mapped BAM file, it looks like you need to use grep on the uncompressed SAM fields:

samtools view -h mapped.bam | grep -v -e 'XA:Z:' -e 'SA:Z:' | samtools view -b > unique_mapped.bam

Explanation follows...

I'm going to assume a situation in which a bioinformatician is presented with a mapped BAM file produced by BWA, and has no way of getting the original reads. One high-effort solution would be to extract the mapped reads from the BAM file and re-map with a different mapper that uses the MAPQ score to indicate multiple mappings.

... but what if that were not possible?

My understanding of BWA's output is that if a read maps perfectly to multiple genomic locations, it will be given a high mapping quality (MAPQ) score for both locations. Many people expect that a read that maps to at least two locations can have (at best) a 50% probability of mapping to one of those locations (i.e. MAPQ = 3). Because BWA doesn't do this, it makes it difficult to filter out multiply-mapped reads from BWA results using the MAPQ filter that works for other aligners; this is likely to be why the current answer on Biostars [samtools view -bq 1] probably won't work.

Here is an example line from a BWA mem alignment that I've just made. These are Illumina reads mapped to a parasite genome that has a lot of repetitive sequence:

ERR063640.7     16      tig00019544     79974   21      21M2I56M1I20M   *       0       0       TATCACATATCATCCGACTCAGCTCGACGAGTACAATGCTAATTTAACACTTAGAATGCCCGGCAATGAAATTCGTTTTCCGTCAATTCTTGAAAATTTC    <AABBEGABFJKKKIM7GHKKJK>JLKLDGMHLIMIHHCGIJKKLJKLNJGLLLKLILKLMFNDLKGHJEKMKKMIJHGLOJLLLKIJLKKJEJLIGG>D    NM:i:4  MD:Z:83A13      AS:i:77 XS:i:67 XA:Z:tig00019544,-78808,21M2I56M1I20M,6;tig00019544,-84624,79M1I20M,6;tig00019544,-79312,33M4I42M1I20M,8;

BWA mem has found that this particular read, ERR063460.7, maps to at least three different locations: tig00019544, tig00019544, and tig00019544. Note that the MAPQ for this read is 21, so even though the read maps to multiple locations, MAPQ can't be used to determine that.

However, the alternative locations are shown by the presence of the XA tag in the custom fields section of the SAM output. Perhaps just filtering on lines that contain the XA tag will be able to exclude multiply-mapped reads. The samtools view man page suggests that -x will filter out a particular tag:

$ samtools view -x XA output.bam | grep '^ERR063640\.7[[:space:]]'
ERR063640.7 16  tig00019544 79974   21  21M2I56M1I20M   *   0   0   TATCACATATCATCCGACTCAGCTCGACGAGTACAATGCTAATTTAACACTTAGAATGCCCGGCAATGAAATTCGTTTTCCGTCAATTCTTGAAAATTTC    <AABBEGABFJKKKIM7GHKKJK>JLKLDGMHLIMIHHCGIJKKLJKLNJGLLLKLILKLMFNDLKGHJEKMKKMIJHGLOJLLLKIJLKKJEJLIGG>D    NM:i:4  MD:Z:83A13  AS:i:77 XS:i:67

... so it filtered out the tag (i.e. the tag no longer exists in the SAM output), but not the read. There are no useful bits in the FLAG field to indicate multiple genomic mappings (which I know can be filtered to exclude the read as well), so I have to resort to other measures.

In this particular case, I can use grep -v on the uncompressed SAM output to exclude alignment lines that have the XA tag (and re-compress to BAM afterwards, just to be tidy):

$ samtools view -h output.bam | grep -v 'XA:Z:' | samtools view -b > output_filtered.bam
$ samtools view output_filtered.bam | grep '^ERR063640\.7[[:space:]]'
<no output>

Hurray! reads filtered. As a little aside, this grep search has a fairly substantial computational load: it's looking for some string with the text XA:Z: somewhere in the line, and doesn't actually capture every situation. Some masochistic person might come along at a later date and decide that they're going to call all their reads HAXXA:Z:AWESOME!:<readNumber>, in which case a tweak to this grep search would be needed to make sure that there's a space (or more specifically a tab character) prior to the XA:Z:.

Now I do a check for any duplicated read names, just to be sure:

$ samtools view output_filtered.bam | awk '{print $1}' | sort | uniq -d
ERR063640.1194
ERR063640.1429
ERR063640.1761
ERR063640.2336
ERR063640.2825
ERR063640.3458
ERR063640.4421
ERR063640.4474
ERR063640.4888
ERR063640.49
ERR063640.4974
ERR063640.5070
ERR063640.5130
ERR063640.5300
ERR063640.5868
ERR063640.6116
ERR063640.6198
ERR063640.6468
ERR063640.6717
ERR063640.6797
ERR063640.7322
ERR063640.750
ERR063640.7570
ERR063640.7900
ERR063640.8115
ERR063640.8405
ERR063640.911
ERR063640.9206
ERR063640.9765
ERR063640.9986

Oh... damn. I wonder what they are:

$ samtools view output_filtered.bam | grep '^ERR063640.3458[[:space:]]'
ERR063640.3458  16  tig00002961 5402    60  58S38M  *   0   0   AGGTACCATTCGATAGAGGGAGAAAGGCACTACTAAAGATTTTGCCACATTTGCTATATCCGTATCGCGAAGATCAGGACTTACTCCGCAGAAGAA    DD6HFFJBKFH=KDILKLGLJEKLKGFJIH8IKHLLMJEK:L:HBGJIHJKFLLKIHJDHLNKCK;KMKGMFKJILIIIMKI9JLKKHEJFII?CC    NM:i:0  MD:Z:38 AS:i:38 XS:i:0  SA:Z:tig00002377,202353,-,14M3I5M1I35M38S,19,5;
ERR063640.3458  2064    tig00002377 202353  19  14M3I5M1I35M38H *   0   0   AGGTACCATTCGATAGAGGGAGAAAGGCACTACTAAAGATTTTGCCACATTTGCTATA  DD6HFFJBKFH=KDILKLGLJEKLKGFJIH8IKHLLMJEK:L:HBGJIHJKFLLKIHJ  NM:i:5  MD:Z:5G48   AS:i:35 XS:i:27 SA:Z:tig00002961,5402,-,58S38M,60,0;

Aha! Supplemental alignments, which use the official SA tag [other canonical alignments in a chimeric alignment]. These appear to be situations where a single read has been split up, and maps to multiple locations. Note that in this case, both of the alignments still have MAPQ scores of over 3. It sounds like the questioner would also want to get rid of these situations as well. This time, there are standard flag fields as well to deal with these situations (0x800: secondary alignment). Except it's not enough to just filter the supplemental alignment, because both read mappings should be removed, rather than just the one (or ones) that happened to be tagged as secondary.

Luckily, BWA appears to put the SA tag into all reads containing supplementary alignments (if this is not the case, I'm sure someone will correct me on that). So, I add in the SA search as an additional grep filter:

$ samtools view -h output.bam | grep -v -e 'XA:Z:' -e 'SA:Z:' | samtools view -b > output_filtered.bam
$ samtools view output_filtered.bam | awk '{print $1}' | sort | uniq -d
<no output>

Done. Easy peasy! </sarcasm>

... that original "high-effort" solution of using a different aligner doesn't look so bad now.


Correlate DEGs from DESeq2, EdgeR and Limma results

I would use the p-value/FDR which each method returns to rank the gene list for that method in order from 'most likely to be DE' to least likely.  You would then have three ranked lists of genes - you would put the actual rankings in for each gene - 1 for the top ranked, 2 for the next, etc.  You could then carry out, say Spearman's rank correlation on the lists of rankings to assess agreement between the methods.

You could also decide to keep all 3 methods, and accept genes as DE if, for example, at least 2 out of 3 methods give an acceptable p-value/FDR.

From my experience, correlation should be reasonably high at the top of the list, where the methods tend to agree more on the rankings, and when the correlation gets lower, further down the list, would be where I would start to question whether the genes are really differentially expressed.


Removing PCR duplicates in RNA-seq Analysis

For normal RNA-seq PCR duplicates are normally kept in, but the duplication rate can be used as a quality control: The higher the duplication rate, the lower the quality. For expression analysis, it is probably best to discard high duplication rate samples, rather than deduplicate them.

In general, the smaller the amount of RNA input into the library prepartion the worst the duplication. Many protocols for very low input quantities (such as single cell) include random barcodes called UMIs (Unique Molecular Identifiers). These allow PCR duplicates to be distinguished from genuinely independent molecules that just happen to have the sample mapping position.


CNSeg no method coercing S4 class to vector

If you are using the latest version of Bioconductor you should always look at the latest documentation of the package in Bioconductor. There could be changes from 2011 to 2018 (In this case there aren't any changes in the section you mention).

The problem is the package doesn't define how dist should behave with that object class (filteredrs is a RS object if I read the source code well).

You can extract the metrics to calculate the distance yourself or ask the maintainer to add this method to the class. However it could be that the maintainer answer you if you post the question in the support site of Bioconductor.


It seems that it is not possible to apply strip() directly on the wildcard placeholders. Luckily it is possible to executed functions in the input directive.

So try changing your input of the complexity_20mer_counter rule to this:

lambda wc: os.path.join(fastq_trimmed_dir, wc.sample + '_' + wc.read + fastq_extension.split('.')[0] + '_val_' + wc.read.strip('R') + '.fq.gz')

The wildcard object is passed to our anonymous function. So have access to the wildcards value via wc.sample and wc.read, which are now handled as normal strings and we can apply strip on them.


How to generate conformations for a ligand for in silico screening?

One option I found was Frog2

https://mobyle.rpbs.univ-paris-diderot.fr/cgi-bin/portal.py#forms::Frog2

https://github.com/tuffery/Frog2

From an input of a mol2 or sdf file it creates a plethora of conformers, which can be outputted as a pdb file


Does Prokka or gff3 change the sequence / subject ID of fasta files?

Prokka is an annotation tool. It takes your contigs file with fasta headers >NODE_108...and searchs for ORFs and their putative functions. Prokka produces submission ready output, the new fasta headers you see >JHIOGLGG_28...are locus tags and gene description from the annotation and are not related to contig node data. Please see the github information page for more info.


Why is Mann-Whitney U-Test over an ExpressionSet yielding different results than original paper

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

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

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

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

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

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

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

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

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

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

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

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

#' Session info
sessionInfo()
```


vcftools: indel size histogram command returns empty file

It looks like lumpy looks for breakpoints in your reference. Presumably places where the is a haplotype switch detectable with the data on-hand?

My only thought is that since lumpy is being used to produce the VCF file that it is inherently not looking for insertions. I would suggest running the same BAM file through the GATK best practices pipeline (easier said than done) and vcftools pipeline. If those produce insertions but lumpy does not, it must be due to how lumpy produces the calls. Good luck!


The problem is that the SRA website can be a little confusing. View the same project on ENA where the organization is more human-readable: https://www.ebi.ac.uk/ena/browser/view/PRJNA574273


Write the output of the VariantAnnotation package to file

Is there any way to force R to give me the same representation as head?

head just cuts off everything after the first six entries. The magic, as it were, happens in the default print method for the specific type. But the representation generated by it is only useful for ad-hoc display, not for storage in a structured data file.

The best solution is therefore to manually convert the column types. Since it’s only one column:

coding$ALT = sapply(coding$ALT, as.character)

I’m a bit puzzled that the REF and ALT column have different types, anyway: that’s probably a bug in the {VariantAnnotation} package that should be reported.

Furthermore, you probably want to use quote = FALSE for the write.table. In general, write.table has terrible defaults. readr::write_tsv is almost always a better choice.


When should I apply a multiple comparison adjustment?

The p-value adjustments corrects the inflated number of significant results when performing many significance tests.
If you have a test with 5 % false positive rate and but you run it several times the chances of getting a false positive results in any of these tests is far greater than 5 %. For example repeating the test with 5 % false positive rate 20 times means that there is a 60 % chance of at least one false positive result.

You can see this in the tutorial you have linked to in your question, where the number of significant differences goes down from 26 without correction to 18 with correction.

Therefore, if you are performing multiple comparisons you should do a correction. The problem of inflated p-values will occur when you are performing many tests even if the tested variables are independent.


Please look up flavivirus 'double loops' as you described them previously (post for "Coronavirus RNA') and associated RNA secondary structure anomalies for dengue virus and associated vaccine (Butantan) and the yellow fever virus and its vaccine (17D). If you are aware of "double loops", we must be aware of the association of RNA secondardy structure and attenuated strains, i.e. live vaccines. The point is that modification of the RNA secondary structure is a key strategy in flavivirus attenuation.


Sequenticon is a Python library for generating identicons for DNA sequences. For instance, the sequence ATGGTGCA gets converted into the following icon:

A web interface is also provided at EGF CUBA (Collection of Useful Biological Apps): Render Sequenticons
Disclaimer: I'm the current maintainer of Sequenticon


I'm not sure that boxplot will be the more appropriate representation, as you will end with two numbers (count of Clonal and SubClonal) per groups of patients.

One solution will be to first create a new categorical variable based on CCF values using for example ifelse statement for example by writing:

df$Clonal <- ifelse(df$CCF > 0.95,"Clone","SubClonal")

Then you can get the count of each Clone and subclone for each responders and non responders by using table:

DF <- as.data.frame(table(df[,c("Response","Clonal")]))

You can finally convert these new values as percentage by grouping according to the Response. For example, using dplyr, you can do something like that:

library(dplyr)

DF %>% group_by(Response) %>% mutate(Freq_percent = Freq / sum(Freq))

And finally, you can plot it as a barchart or single point in ggplot2.
Hope it helps you to figure it out how to deal with your data.


How to calculate logCPM across all samples?

If you want the values by group then subset y to contain the samples of interest and then feed that to aveLogCPM().


Why ORF13 and ORF14 of Bat Sars coronavirus Rp3 have no corrispondence in Sars2?

The virus Rp3 stands for Rhi pea isolated in China in 2004.
The trees is here and groups with the Bulgarian sarbecoviruses and is closely related to the African sarbecoviruses. Thus although it is a Chinese sarbecovirus is not with the "Chinease originating" clade as is SARS and Covid-19.
Thats why the more mysterious ORFs ORF13 and 14 are missing, because over notable genetic distances these genes along with ORF10 are present/absent. This is a common feature of the betacoronaviruses.
In summmary this virus is notably distance from the major human pathogenic SARS/covid-19 viruses. Thus the absence of these auxilary genes isn't really surprising.
Its the very last taxa on the SARS-related viruses clade (monophyly) and is in a different group to SARS, which is the Beta-1 clade on this most excellent paper.

To answer the comment. They are outgroups, basically ancestral to COVID-19 and SARS-1 clade. So they lost these genes whilst the ingroup including COVID-109 and SARS-1 retained ORF13, 14 ... but I think lost ORF10 from memory


I have seen many posts regarding counts to RPKM and TPM.

There’s your answer then: FPKM = RPKM. It’s simply a more accurate name.

Speaking of RPKM for paired-end data is discouraged because the reference to “read” in this context lends itself to ambiguity. But mathematically the quantity is the same: we are counting fragments, not individual reads (of which each fragment has two, for paired-end data).

But as mentioned in the comments, there are good reasons against using FPKMs.


If Python is an option, you could check BioPython and GenomDiagram. This should get the job done.


How to calculate the number of peaks that are upstream/downstream of some other peaks

You can use windowBed from bedtools with -c if you just want counts. Or use -wa -wb and then use groupBy for any other operations.


Question: Bismark methylation for non-uniquely mapping BS-seq reads

I'm not certain about your piRNA approach, maybe others will be able to help with that question.
As for multi-mapped reads, from Bismark's documentation you can find the --ambig-bam flag:

For reads that have multiple alignments a random alignment is written out to a special file ending in .ambiguous.bam. The alignments are in Bowtie2 format and do not any contain Bismark specific entries such as the methylation call etc. These ambiguous BAM files are intended to be used as coverage estimators for variant callers.

It doesn't look like Bismark will return all of the multiple alignments, just a randomly selected one for a given read or read pair (this probably results from the way Bowtie2 handles multiple alignments within Bismark).


How can I use Arlequin via the command line?

For MAC OSX, the executable found in arlsumstat_macosx (arlsumstatmac_64bit) seem to work appropriately. The directory for arlsumstat_macosx does not contain any examples but you can use example files from arlecore_macosx/Example files_linux using this syntax.

The complication is that arlsumstatmac_64bit will look in the current directory for the file arl_run.ars. There is a file with the same name in arlecore_macosx but it is not compatible. So you have to be in the arlsumstat_macosx.

In short, cd to arlsumstat_macosx and do

./arlsumstatmac_64bit path/to/arlecore_macosx/Example\ files_linux/Freqncy/cohen.ars path/to/arlecore_macosx/Example\ files_linux/Freqncy/cohen.arp

, where the .arp file contains the data and the .ars file contains the settings (the description of what you want Arlequin to do).


How to set new identities for new clusters based on a gene expression in Seurat?

You have to define your ident it has to have the structure of a table. You put the cell names as names and clusters as value (for instance cell 1 = YEPF+ cell 2 = YEFP-):

pident=as.factor(clusters)
names(pident)=cellNames
object1@ident=pident

Then you can run the FindAllMarkers function on this object:

FindAllMarkers(object1, min.pct = 0.25, min.diff.pct = 0.25)


Tools or softwares for clustering evaluation

The definitive answer is no, that is why they are not so popular. There are metrics that can be used:

you can examine "minimum distance" as a criterion under the assumption that the lowest distance between taxa represents the largest homology (e.g. if it is sequence data), which would be translate to other biological scenarios.
resampling, e.g. bootstrapping, jack-knifing and examine the tree structure of each node against the biological plausibility of the answer being correct.

For both points a biological a priori, i.e. biological plausibility of assessing the tree structure is required to confirm/reject the clustering method.

Clustering algorithms are not like the probabilistic methods, where the highest likelihood is the singular criteria and of course they do give different answers. Its a brief explanation I know but hope it makes sense.


merge rows (to columns) in R dataframe based on IDs

For a data.frame named foo:

sink('some_output_file.txt')
uvals = unique(foo[,1])
for(uIDX in c(1:length(uvals))) {
    gname = uvals[uIDX]
    IDX = which(foo[,1] == gname)
    cat(sprintf("%s\t%s\n", gname, paste0(foo[IDX,2], collapse="\t")))
}
sink()


Convert paired-end BAM into a single-end BAM and keep all the reads

Like Devon pointed out, most likely you should sort out whether the files have been marked for duplicates correctly.

You can also use samtools rmdup too

samtools rmdup <bamfile> <output>

Also, I find it very odd that 50% of the reads are lost when you keep only properly paired reads.

I would suggest removing the duplicates and then converting the bam file to a bed file. Macs2 can take in the bed format. Any extremely crude python script that converts a bam into a bed:

import pysam
samfile = pysam.AlignmentFile("test.bam", "rb")
open("test.bed", 'w').close()
f = open("test.bed", "a")
STRAND = ["+","-"]
for read in samfile.fetch():
    STR = STRAND[int(read.is_reverse)]
    BED =[read.reference_name,read.pos,read.reference_end,".",read.mapq,STR]
    f.write('\t'.join([str(i) for i in BED]))
    f.write('\n')
f.close()

You can also try the conversion with bedtools or definitely refine that script above with pandas etc.


ethnicity check either from bam or vcf files

The main difficulty here is the use of GRCh38.  Unfortunately, despite the fact that it's more than four years old, the major ethnicity-labeled public datasets (1000 Genomes, gnomAD when allele frequencies are enough) still aren't available for that reference.  It is necessary to perform a liftover operation, or just use overlapping rsIDs and hope for the best.

Let's suppose you go with overlapping rsIDs, and a list of those rsIDs, one per line is in 'rsids.txt'.  (Due to strand flips between reference builds, you may want to restrict this list so that all allele codes match, and it's also reasonable to throw out A/T and C/G SNPs.)  Then, the following process would work:

Download ADMIXTURE (https://www.genetics.ucla.edu/software/admixture/download.html ), plink (https://www.cog-genomics.org/plink/1.9/ ), and/or chr1-chr22 of 1000 Genomes phase 3 (ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ ) if you don't already have them.
Use plink to extract the overlapping higher-MAF rsIDs from each 1000 Genomes phase 3 VCF ("plink --vcf ... --extract rsids.txt --maf 0.05 --make-bed --out ..."), and then merge the resulting per-chromosome filesets ("plink --merge-list ... --out merged_phase3_subset").
ADMIXTURE prefers a dataset with around 100k variants.  Unless you have a small list of overlapping rsIDs (in which case you had better use another method), you should have more than that.  plink's LD-pruning function is a good way to select a subset for ADMIXTURE's use:
"plink --bfile merged_phase3_subset --indep-pairwise 500kb 1 0.2; plink --bfile merged_phase3_subset --extract plink.prune.in --make-bed --out admixture_data"

Adjust the 0.2 threshold as necessary to keep the right number of variants.

Run ADMIXTURE in unsupervised mode ("admixture admixture_data.bed 5 -j8"; adjust the -j parameter depending on the number of processor cores).  This generates an admixture_data.5.P file with population allele frequencies, and an admixture_data.5.Q file with sample population assignments.  Verify that the .Q file corresponds to the 1000 Genomes phase 3 'superpopulations'.
Convert your data to plink-format if necessary, keeping just the overlapping rsIDs and 'downgrading' to GRCh37 coordinates.  Make sure your variants are sorted in GRCh37-coordinate order, and the allele order also matches that in your reference dataset ("--a2-allele admixture_data.bim 6 2" during your final --make-bed operation will do the trick); otherwise ADMIXTURE won't do the right thing.  Then run ADMIXTURE in projection mode ("cp admixture_data.5.P my_data.5.P.in; admixture -P my_data.bed 5").  my_data.5.Q will then have the ethnicity estimates you're looking for.


Can we use non base-called fast5 files in poretools?

No, poretools does not do basecalling. The poretools fastq command can be used to extract the FASTQ information from the basecalled FAST5 file (via MinKNOW live-basecalling or albacore). Alternatively, both of these basecallers can export a FASTQ file directly if desired.

The difference between a basecalled FAST5 and a nonbasecalled FAST5 is the presence of the FASTQ within the file. This can be directly viewed using an HDF5 viewer, like h5dump from the commandline, or HDFView from the GUI.


Error When Using biocLite as an installer in rpy2 python library

You're using an outdated version of R. You can use the R 4.0 from conda-forge channel, just create an environment with r-base and then install rpy2, the latest version compatible with R 4.0 is available through pip
conda create -n rpy2 -c conda-forge r-base
conda activate rpy2
pip install rpy2

Now you can use rpy2 with R 4.0 in your python code.
import rpy2.robjects as robjects
robjects.r('install.packages("BiocManager", '
           'repos="http://cran.r-project.org")')
robjects.r('BiocManager::install("seqLogo")')
robjects.r('library(seqLogo)')

Remember deactivate your env with conda deactivate.


After running a DRIMSeq pipeline, how do I know which genes are upregulated in the different conditions?

After running the test for differential transcript usage you get a list with gene- and transcript-level p-values that you can use to filter which differences are significant. The gene-level p-value tells you if there is differential usage of transcripts within a gene and the transcript-level p-value tells you wich transcripts actually are changing.  Therefor you will want to select significant genes for further investigation.

The plots are mainly used to visualize top hits as detailed in the vignette.

res <- results(d)
# The data is ordered by p-value
res <- res[order(res$pvalue, decreasing = FALSE), ]
# The top hit is selected
top_gene_id <- res$gene_id[1]
plotProportions(d, gene_id = top_gene_id, group_variable = "group")

However you should consider running the two-stage test that is described on page 16 of the vignette you are using to reduce the number of false positive hits. From this you will get another list of adjusted p-values that can be used as above.


awk find Sickle cell anemia from father, mother, son, daughter


Continuing analysis based on a log-likelihood ratio? (Comparing methylation between two PoolSeq phenotypes)

I assume you don't want to spend 1-2 years developing your own method based on nanopore data, so the path forward is to convert the output into a format used by standard DMR calling programs. There's no standard format for that, annoyingly, but the general things you'll need are:

Coordinates of a C/CpG with coverage
Number of reads supporting methylation at that site.
Number of reads not supporting methylation at that site.

Typically you would do some filtering before/after making such a table so you only have decent quality calls (that would translate into some reasonable log likelihood for your data) and possibly only sites with decent (5x, 10x, 30x, etc.) coverage. You'd then feed that for all of the replicates in all of your groups into programs like bsseq, DSS or Metilene (to not so randomly select 3 examples).

If there aren't DMR finding programs already developed specifically for nanopore data I expect someone is developing it. The technology has a higher error rate, but given that bisulfite conversion as required in Illumina data also leads to a higher error rate that's not necessarily a huge issue if the methylation calling itself is reasonably accurate and the depth is high enough. Note that you probably need 10X coverage as a bare minimum unless your effect size is very very large.


What does "fetching by region is not available for SAM files" mean?

That isn't actually a bam file as John Marshall figured out. I am keeping the rest of my answer since it could be useful to someone else, but the issue here was that you had a compressed (bgzipped) sam file and not an actual bam file and that's why you were getting that error. When I sorted your file in preparation for indexing it, I converted to a bam which is why the rest of my answer worked.

You don't have the index file for the bam file you're using. I used this script on the file you linked to (changing the name so that it corresponds to the right file and a contig in that file):

#!/usr/bin/env python3

import pysam
samfile = pysam.AlignmentFile("3_Tce_1_mapped.bam", "rb")
for read in samfile.fetch('3_Tce_b3v08_scaf005149'):
    print(read)
samfile.close()

The directory I ran it in had:

$ ls 3*
3_Tce_1_mapped.bam

And I got the error you described:

$ foo.py
Traceback (most recent call last):
  File "/home/terdon/scripts/foo.py", line 5, in <module>
    for read in samfile.fetch('3_Tce_b3v08_scaf005149'):
  File "pysam/libcalignmentfile.pyx", line 1107, in pysam.libcalignmentfile.AlignmentFile.fetch
ValueError: fetching by region is not available for SAM files

However, indexing the bam file fixed it:

$ samtools sort 3_Tce_1_mapped.bam > 3_Tce_1_mapped.sorted.bam
$ mv 3_Tce_1_mapped.sorted.bam 3_Tce_1_mapped.bam
$ samtools index 3_Tce_1_mapped.bam
$ ls 3*
3_Tce_1_mapped.bam  3_Tce_1_mapped.bam.bai

$ foo.py | wc
227    2724   16725

So just sort and index your files before attempting to seek in them. Which makes sense since the index's job is primarily to allow seeking.


Converting gene names from one public database format to another

I tend to use Ensembl Biomart for such queries since there are APIs for various programming languages, e.g. biomaRt, and, maybe more interestingly, via a REST API (although it’s a pretty terrible one).

To translate identifiers from different databases, proceed as follows:

Choose database “Ensembl genes”
Choose dataset your desired oganism
Go on “Filters” › “Gene:” › “Input external reference ID list”

Select the chosen source database
Provide a list of IDs, delimited by newline

Go to “Attributes” › “Gene:” › untick “Transcript stable ID”

If Ensembl IDs are desired, leave “Gene stable ID” ticked …
Otherwise untick it; go to “External:”, tick your desired identifier format

Click “Results” at the top left. This gives a preview that can be exported into various formats; alternatively the top centre buttons “XML” and “Perl” provide the query in XML (for SOAP/REST requests) and as a (horrendously formatted) executable Perl script.


General framework for fusing biological networks

I believe you're referring to something related to multi-layer networks. Here you can find an interesting reading about this kind of networks:

https://academic.oup.com/comnet/article/2/3/203/2841130

Also you can try with this one:

https://www.sciencedirect.com/science/article/pii/S0370157314002105

As a matter of fact, I'm not an expert of multi-layer networks but, as far as I can tell, the main problem (as you mentioned) is to find correlations, or connections, between the different layers in order to be able to draw some meaningful conclusion.

I am not sure you can call this a multi-layer network but this can be an example of what you're looking for: suppose you have a protein-protein interaction network. You can enrich your network using biological information. You can associate biological functions, or cellular locations or whatever you can find in GeneOntology (just to make an example of a biological database), as a further layer of information regarding your protein-protein interaction network.

This will help you in elucidating the role a node (or more nodes) is playing with respect to its neighbours: let's say five nodes form a clique and you also know they share the same cellular location...then you're able to say something about these nodes (maybe they form a complex?) that has a dual nature: biological and graph-based.

Hope it helps!


HISAT2 - are all multi mapped reads reported by default?

You said:

  ...the first quote means that HISAT2 finds all distinct valid alignment and reports them by default.

Not exactly. The first section you quoted doesn't say it will find and report all distinct alignments, but rather that it would find and report one or more distinct alignments. It's a subtle difference, but an important one.

That section describes how distinct mappings are defined, not how many are reported. The second section you quoted, as well as the documentation for the -k parameter, are all consistent with this interpretation.

-k <int>        It searches for at most <int> distinct, primary alignments for
                each read. Primary alignments mean alignments whose alignment
                score is equal or higher than any other alignments. The search
                terminates when it can't find more distinct valid alignments, or
                when it finds <int>, whichever happens first. The alignment
                score for a paired-end alignment equals the sum of the alignment
                scores of the individual mates. Each reported read or pair
                alignment beyond the first has the SAM 'secondary' bit (which
                equals 256) set in its FLAGS field. For reads that have more
                than <int> distinct, valid alignments, hisat2 does not guarantee
                that the <int> alignments reported are the best possible in
                terms of alignment score. Default: 5 (HFM) or 10 (HGFM)

                Note: HISAT2 is not designed with large values for -k in mind,
                and when aligning reads to long, repetitive genomes large -k can
                be very, very slow.

By default, 5 or 10 distinct mappings are reported (depending on the mode), but you can control this with the -k parameter.


Does SBOL require representing intermediate products like mRNA?

As far as I understand, the purpose of SBOL is not to create a comprehensive model, but to communicate a design.  So the question you have to ask is "is the mRNA part of my design, or is it an implicit part of the biology which does not help to convey what I'm trying to build?".

I think your summary is a giveaway:

  For example, let's say that I'm representing transcriptional regulation of the pTet promoter by means of a coding sequence that expresses TetR.

Notice that the mRNA was not mentioned?  I also daresay that if you drew this design as a diagram, the mRNA would be omitted - another indication that it is not helpful in communicating what you designed.


Is it possible to go from a digital sequence to an actual nucleotide?

Oligonucleotide synthesis is a service offered by many providers such as IDT.


How to make a distinction between the "classical" de Bruijn graph and the one described in NGS papers?

Several papers have made this distinction, and a few indeed use different terms to distinguish between them. For example, Kazaux et al. (2016) acknowledge that:

  These constraints favour the use of a version of the de Bruijn Graph (dBG) dedicated to genome assembly – a version which differs from the combinatorial structure invented by N.G. de Bruijn.

Kingsford et al. (2010) also recognise the distinction:

  Note that this definition of a de Bruijn graph differs from the traditional definition described in the mathematical literature in the 1940s that requires the graph to contain all length-k strings that can be formed from an alphabet (rather than just those strings present in the genome).

The oldest reference I found for a specific term to refer to the assembly-related structure is Skiena and Sundaram (1995), where they call it a subgraph of the de Bruijn digraph. Later, in 2002, Błażewicz et al. will refer to it as a de Bruijn induced subgraph.  The term de Bruijn subgraph is also formally defined in Quitzau’s thesis (2009). There, and also in the article (Quitzau and Stoye, 2008) the authors describe the sequence graph as a modification of the sparse de Bruijn subgraph (commonly used in assembly problems), where non-branching paths are replaced by a single vertex. The term sparse de Bruijn graph is also used by Chauve et al. (2013).

Another term that I found was word graph, described by both Malde et al. (2005) and by Heath and Pati (2007) as a subgraph or as a generalization of a de Bruijn graph. Rødland (2013) summarises some of the terms used for this data structure:

  The data structure is best understood in terms of the de Bruijn subgraph representation of S[k]. (...) Some authors may refer to this as a word graph, or even just a de Bruijn graph.

Although we can recognise that the distinction is not very relevant, the question is asking specifically for the situation where one wants to make such a distinction.


How do you query and explore ENCODE data?

I don't know whether there is an API, but ENCODE's website does provide an interactive data matrix where you can filter data based on assay and sample type, place data sets in a "shopping cart", and then proceed to "checkout" to download the files of interest.


How can I calculate loss of heterozygosity (LOH) in NGS sequencing data?

Answer from @fabio-marroni, converted from comment:
WARNING! Self-promotion! We published a freely-available tool for detection of LOH (and also redeuction fo homozygosity). Here you find the link to the paper and therein will find the link to the software. We compared with other tools that you might also try to see which one is the best performing. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5466433/


Interpreting a kernel density function plot of gene expression

There are two topics here to be discussed. Gene Set Enrichment Analysis and a kernel Density plot.
GSEA
The Wikipedia page about Gene Set Enrichment Analysis covers well the principle and the tools that can be used to calculate it. To cite it:

The general steps are summarized below:
Calculate the enrichment score (ES) that represents the amount to which the genes in the set are over-represented at either the top or bottom of the list. This score is a Kolmogorov–Smirnov-like statistic [a ranking based statistic].
Estimate the statistical significance of the ES. This calculation is done by a phenotypic-based permutation test in order to produce a null distribution for the ES. The P value is determined by comparison to the null distribution.
Calculating significance this way tests for the dependence of the gene set on the diagnostic/phenotypic labels

Basically, you want to test if the ranking is random or not. Therefore you make a random dataset and test whether your set is like that. One assumption underlying this is that the output of the gene whatever it may be is not normally distributed (parametric, follows a bell curve), hence the use of ranking, which removes this (you lose signal but you lose biases). I am assuming the non-parametric nature of your ouput values was tested or is a known non-parametric value, such as RPKM or log2 fold change (by virtue of being a ratio). In other words, the magnitude difference between your samples is lost, but their order is not.
There are several tools, some online, that do this enrichment testing for you depending on the task, some even work with gene ontologies.
KDE
A kernel density plot is a like a histogram, but smoothed, albeit not in a moving average way. In a histogram you divide your x-axis into bins: it is discreet and you get an integer count per bin. In a kernel density plot the data is fitted to "probability density function", an equation which given x will give y, where the integral of the curve is one, hence why y is density not counts. The x axis will be your bins-like value, in this case rank. The binning window metric for a KS function is called bandwidth.
See Matlab documentation for a fairly non-opaque explanation.
However, if you want to compare histograms, I strongly suggest violin plots, not because they are strongly in vogue, but because you can compare histograms without worrying about overlaying them.


ribo interval list for picard ensembl mm9 gtf

At least with Ensembl annotations, you can grep rRNA foo.gtf > rRNA.gtf to get the annotated rRNAs. Note however, that the 45S cassette isn't in the mouse reference genome, so you'll never have alignments to 18S, 5.8S, or 28S rRNA. I always get the feeling that the Collect*Metrics programs are made with human sequences in mind, since the cassette is included in the human reference genome.

Anyway, while you can use the annotation and reference sequence to make picard happy, make sure to append the Rn45S sequence to your genome if you want to actually measure rRNA levels. In my experience this is only useful for rRNA-depleted samples, since the rates can vary wildly there.

Edit: To make a file compatible with picard's "interval list", one can do the following with the GTF file:

grep rRNA foo.gtf | awk 'BEGIN{OFS="\t"}{print $1, $4-1, $5}' > foo.bed

According to the picard documentation a BED file can be used as input and the above awk code extracts the appropriate columns and switches to 0-based coordinates used in the BED format.


You can omit --nodes, you need the following:

#!/bin/bash
#SBATCH --ntasks-per-node=1
#SBATCH -c threads_from_snakemake
#SBATCH -p some_partition

Commands go here

For slurm, you may want to modify my SlurmEasy script. We use that with snakemake all the time, usually of the form:

snakemake -c "SlurmEasy -t {threads} -n {rule} --mem-per-cpu {cluster.memory}" ...

You'll also want units with your memory requests, like 12G. Also, in general it's best not to submit snakemake itself to the cluster, rather run it on an interactive node and have it submit jobs on your behalf (if the admins complain point out that snakemake is barely using any resources on the head node).


The kind of normalization depends on what you what to explore. So, there is no absolute answer to this. Different normalizations highlight different aspects of your dataset. There is a nice paper for a higher throughput quantitative study that I would recommend reading:
Nusinow, David P., John Szpyt, Mahmoud Ghandi, Christopher M. Rose, E. Robert McDonald 3rd, Marian Kalocsay, Judit Jané-Valbuena, et al. 2020. “Quantitative Proteomics of the Cancer Cell Line Encyclopedia.” Cell 180 (2): 387–402.e16.
What you want to aim for is a normalization that removes batch effects and that clusters together your quality control samples that you know are very similar to each other.


Extract sequences from MultiFASTA aligned file, by coordinates

You can use SeqIO for the writing in addition to reading the sequences. This has the advantage that it will always write the correct fasta header, which is not always just the .id if there is whitespace, e.g. (untested)
from Bio import SeqIO

with open("Spike protein.fasta", "w") as out:
    for seq_record in SeqIO.parse("Whole genome.fas", "fasta"):
        seq_record.seq = seq_record.seq[21563:25384]
        SeqIO.write(seq_record, out, "fasta")

BTW don't put spaces in your filenames, it can lead to issues


wtdbg2: practical implications of k-mer fsize and psize choice

There is some description of homocompressed polymers (-p kmer) in this paper: https://arxiv.org/pdf/1708.01492.pdf


In silico PCR: any models that can quantitatively predict the initial content from the end amplified result?

The basic equation for PCR kinetics (Equation 1) states that the amount of amplicon after c cycles (Nc) is the starting concentration of the amplicon (N0) times the amplification efficiency (E) to the power c. The PCR efficiency in this equation is a number between 1 and 2 (2 indicates 100% efficiency).  NCBI

In this  article  on research gate they said  amplifx  was the software you were looking for to, calculate "(i.e. expected efficiency in PCR reaction)". The amplifx even mentions "Testing primers in in silico PCR".


How do I get repeating regions from fasta file or directly from internet

Click on the 'GenBank' dropdown menu at the top left corner of this page to see a graphical view of the sequence. You should see a 'Repeat region' track in the Sequence Viewer. Click on the 'Download' button in the top right corner of the Sequence Viewer, choose 'Download Track Data' to launch a dialog box where you can enter a range, select the 'Repeat region' track and download data in BED format. Unfortunately, the names for the features are not great.
If you want the names to be the same ones in the flat file, you can download the feature table and parse it. To download the feature table, click on the 'Send To' link at the top right corner of this page, choose 'Complete Record', 'File' and 'Feature Table' as 'Format'. After you download the feature table, you can parse it using awk as shown below to create a 4-column table with the data.
cat seq.ft \
  | awk 'BEGIN{FS="\t";OFS="\t"}{
    if(NF==3) {print s,e,f,q ; s=$1; e=$2; f=$3; q=""}
    else if(NF==2) {print s,e,f,q ; s=$1; e=$2; q=""}
    else if(NF==5) {q=q"|"$4"="$5} }' \
  > feat_table.tsv


Is there any way to align ChIP-seq reads to telomeres?

This paper https://www.biorxiv.org/content/10.1101/728519v1.full examines the mapping of some human telomeric sequences, and indicates that there is something like 130kb missing at the telomeres in the HG38 assembly.  One option might be to take the HG38 assembly, and add some missing telomere sequences as an additional pseudochromosome to the FASTA file.  Then if you used this modified reference to align your reads, you should soon get a feel for the level of telomeric reads in your data set.


Pathifier bioconductor package recommends to use technical noise for an argument. How to determine it?

Answer from @llrs, converted from comment:
Probably is a question for the maintainer but I would guess that is a variance of the gene in the control dataset. So, if in normal samples the gene has a standard deviation of 0.5 this is the expected too in the tumoral cells. calculate the sd only on the control samples. So if the healthy/control samples have such variation the altered ones should have these or more according to the article (if I understand it correctly).
std <- apply(data[, controls], 1, sd, na.rm = TRUE); quantile(std, 0.25) should be much faster and efficient
But  I would check it in support.bioconductor with the authors/maintainers.


you can do samtools view pathtobamfile, this will convert the file to .sam format and send it to stdout, the columns can then be parsed out which contain the information you are after.


Difference between paired-end, mate-pair and long read

They are all very different in separate regards, but they all refer to different wet-lab and sequencing protocols/technologies.
First, PE (paired end) reads are typically short (50-300) reads, most often Illumina HiSeq, MiSeq or NovaSeq protocols. Both pairs originate from a single fragment which is sequenced from either end:
------> PE1
-------------------------------- Fragment
                         <------ PE2
       ^^^^^^^^^^^^^^^^^^        Insert

In contrast mate pairs arise from a fragment that is circularized before sequencing. Now imagine the fragment is circularized:
                     -----> MP1
....--------------------------------.......
          MP2 <-----
                    ^Circularized (joined) here

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


How to find novel transcripts using GFFcompare?

The output of gffcompare includes several files per run (just like cuffcompare). Example for a run:

$ ls cuffcmp.* | sed 's/\t/\n/'
cuffcmp.combined.gtf
cuffcmp.loci
cuffcmp.output.gtf.refmap
cuffcmp.output.gtf.tmap
cuffcmp.stats
cuffcmp.tracking

From my experience, the easiest file to manipulate is the tmap file:

  Tab delimited file lists the most closely matching reference transcript for each query transcript. There is one row per query transcript, and the columns are as follows:

1     Reference gene name    The gene_name attribute of the reference GTF record for this transcript, if present. Otherwise gene_id is used.
2     Reference transcript id     The transcript_id attribute of the reference GTF record for this transcript
3     Class code  The type of relationship between the query transcripts in column 4 and the reference transcript (as described in the Class Codes section below)
4     Query gene id   The query (e.g., Stringtie) internal gene id
5     Query transcript id     The query internal transcript id
6     Number of exons     The number of exons in the query transcript
7     FPKM    The expression of this transcript expressed in FPKM
8     TPM     the estimated TPM for the transcript, if found in the query input file
9     Coverage    The estimated average depth of read coverage across the transcript.
10    Length  The length of the transcript
11    Major isoform ID    The query ID of the gene's major isoform
12    Reference match length  The length of the longest overlap with a reference, ‘-’ if there is no such exonic overlap

The important column to look at for transcript classification is actually the 3rd, "class code". The content is meaningful if you used option -r when running gffcompare.

In the tmap file you can access info on both transcript and gene ID directly, both from the assembler and the reference annotation.

Now what you need is:

n. assembled genes
n. novel genes
n. novel transcripts
n. transcripts matching annotation

At this point, it is essential to understand what the authors mean by "novel" (gene/transcript), does it mean "non previously annotated" (= class code u, for "unknown")? Also, what about "matching", does it mean "perfectly matching" (= class code =, for "complete, exact match of intron chain")?

Under these assumptions, you can do some easy stats with some bash help.

n. assembled genes:

$cat cuffcmp.output.gtf.tmap | sed "1d" | cut -f4 | sort | uniq | wc -l
n. novel genes:

$ cat cuffcmp.output.gtf.tmap | awk '$3=="u"{print $0}' | cut -f4 | sort | uniq | wc -l
n. novel transcripts:

$ cat cuffcmp.output.gtf.tmap | awk '$3=="u"{print $0}' | cut -f5 | sort | uniq | wc -l
n. transcripts matching annotation:

$ cat cuffcmp.output.gtf.tmap | awk '$3=="="{print $0}' | cut -f5 | sort | uniq | wc -l

I hope this helps. I think it's an easier approach than trying to "translate the numbers related to novel exons and novel introns to the number and identities of the novel transcripts/genes" from the stats file.

Although you have to think a bit what "class codes" to include in each count. Maybe my assumptions are too strict for your use case.


Invest some time in learning unix tools such as grep, sed and awk.

sed -r '/^>/s/:[0-9]+-[0-9]+//' input.fasta > output.fasta

This line of code matches the pattern :<numbers>-<numbers> on lines that start with > (FASTA sequence headers) and replaces that pattern with '', essentially removing text that matches the pattern.


Biology behind PCA analysis based on SNP

First of all, PCA is a technique for dimension reduction. Basically, the goal is to compare tens of thousands of SNPs in Drosophila. Now if you only have 2 SNPs, you can plot them on a 2D scatter plot. If you have 3 SNPs, you may try a 3D plot. But now imagine you have 30,000 SNPs, but you CANNOT plot a 30000-dimensional plot. To visualize this high dimensional data, what we can do is to perform dimensional reduction like PCA. PCA tries to find a set of orthogonal coordinations that explains most of the variation in the data (if there no variation, there is no information contained in the data, which essentially means there is no data). The idea is that PC1 carries most variation can be explained, and PC2 carries the second most... For lower PCs like PC50 or PC60, they probably only carry noise in the data. Therefore, the higher PCs (PC1, PC2 and so on) effectively summarizes the useful information in the data. So you can visualize the "structure" of the data in a 2D PCA plot.
By looking at the distance between points on a PCA plot, you can tell how similar the two data points are. But if you see two populations that are perfectly separated on PCA plot, it does not mean that the 2 population differ completely at every SNP, because PCA is a summarization of all SNP included.


I was still puzzled from the answers, so I tried to calculate with all the steps. I take this definition "$C_k$ is the number of reads containing a k-mer." and corresponding definition for coverage ($C$): "$C$ is the number of reads covering a base".

Coverage is $C = \frac{T \cdot R}{L}$, where $T$ is total number of reads, $R$ is read length and $L$ is length of genome. Given the $C_k$ definition, $C_k = \frac{T (R - K + 1)}{L-K+1}$, where $R - K + 1$ is just number of kmers in a read, and $L-K+1$ is number of kmers in a genome. Then,

$$C_k = \frac{T (R - K + 1)}{L-K+1} = \frac{T (R - K + 1)}{L-K+1} \cdot \frac{R}{R} = \frac{R - K + 1}{R} \cdot \frac{T \cdot R}{L - K + 1}$$

since $L >> K$, we can approximate $L - K + 1 \approx L$, then we reduce the expression to

$$\frac{R - K + 1}{R} \cdot \frac{T \cdot R}{L} = \frac{R - K + 1}{R} \cdot C$$

which is the formula for $C_k$.