How to make chromosome color maps for bed ranges

You can use karyoploteR to plot your regions on an ideogram quite easily.

Disclosure: I'm one of the authors of the tool.

For this example I'll start creating 4 GRanges with random regions, but with your data you should load your bed into R and split it into four GRanges.

library(karyoploteR)

typeA <- createRandomRegions(nregions = 100, length.mean = 5e6, length.sd = 3e6)
typeB <- createRandomRegions(nregions = 20, length.mean = 15e6, length.sd = 10e6)
typeC <- createRandomRegions(nregions = 1000, length.mean = 5e5, length.sd = 2e5)
typeD <- createRandomRegions(nregions = 100, length.mean = 2e6, length.sd = 5e6)

colors <- c(A="#FFD700", B="#00BFFF", C="#FF7F50", D="#90EE90")

And we can plot them. We have to start creating an empty ideogram with plotKaryotype (we are working with human here, but you can use any genome you want) and then plot the regions in the GRanges with kpPlotRegions. As a bonus, if two regions overlap it will plot them one above the other.

kp <- plotKaryotype(genome = "hg19")

kpPlotRegions(kp, data=typeA, r0=0, r1=0.2, col=colors["A"])
kpPlotRegions(kp, data=typeB, r0=0.25, r1=0.45, col=colors["B"])
kpPlotRegions(kp, data=typeC, r0=0.5, r1=0.7, col=colors["C"])
kpPlotRegions(kp, data=typeD, r0=0.75, r1=0.95, col=colors["D"])
legend(x = "right", legend = names(colors), fill = colors, border=NA)

And you'll get something similar to this

In this example we have fairly large regions and not a lot of them. As @heathobrien pointed out, plotting the actual regions in the whole genome plot might not be the best visualization option. You can either use zoom to create different plots of smaller parts of the genome or change your visualization to show the density of regions along the genome with kpPlotDensity.

Plotting the density and tweaking a bit the visualization (and changing the example to have 10 times more regions) with this codei

typeA <- createRandomRegions(nregions = 1000, length.mean = 5e6, length.sd = 3e6)
typeB <- createRandomRegions(nregions = 200, length.mean = 15e6, length.sd = 10e6)
typeC <- createRandomRegions(nregions = 10000, length.mean = 5e5, length.sd = 2e5)
typeD <- createRandomRegions(nregions = 1000, length.mean = 2e6, length.sd = 5e6)

kp <- plotKaryotype(genome = "hg19", plot.type = 4,
                    ideogram.plotter = NULL, labels.plotter = NULL, main="Density")
kpAddCytobandsAsLine(kp)
kpAddChromosomeNames(kp, srt=45)

kpPlotDensity(kp, data=typeA, r0=0, r1=0.2, col=colors["A"], window.size     = 10e6)
kpAddLabels(kp, labels = "Type A", r0=0, r1=0.2)
kpPlotDensity(kp, data=typeB, r0=0.25, r1=0.45, col=colors["B"], window.size = 10e6)
kpAddLabels(kp, labels = "Type B", r0=0.25, r1=0.45)
kpPlotDensity(kp, data=typeC, r0=0.5, r1=0.7, col=colors["C"], window.size = 10e6)
kpAddLabels(kp, labels = "Type C", r0=0.5, r1=0.7)
kpPlotDensity(kp, data=typeD, r0=0.75, r1=0.95, col=colors["D"], window.size = 10e6)
kpAddLabels(kp, labels = "Type D", r0=0.75, r1=0.95)

it would look something like this

The plots could be tweaked a lot changing various parameters and you could add other data on top of it if needed. You can find more information on the Bioconductor landing page or at the karyoploteR's Tutorial and Examples page.


Transform and feed data into DESeq2 with DESeqDataSetFromMatrix

Error in DESeqDataSet(se, design = design, ignoreRank) : all variables
  in design formula must be columns in colData

I think that means just what it says...you don't have a "Groups" column in coldata.  Did you spotcheck your ColData data frame to make sure it looks the way it should?

> cells <- c("1","2","3")
> groups <- c("a", "b", "c")
> test <- data.frame(Cells <- cells, Groups <- groups)
> test
  Cells....cells Groups....groups
1              1                a
2              2                b
3              3                c
> colnames(test)
[1] "Cells....cells"   "Groups....groups"


Python API for working with ENSEMBL genomes

Now there's another option, the ensembl_rest module, a thin wrapper around the Ensembl REST API to simplify its usage and make it more pythonic. You can find the documentation here.

To clarify things, I'm the creator and maintainer, but still think it's a legitimate alternative.


Calculating polarity and hydrophobicity of a ligand/small molecule

You can calculate their polar surface area or logo using rdkit or openbabel. Other molecular properties may be available as well.


Using Entrez.efetch() to retrive .fasta file from any NCBI database?

In order to download the FASTA sequence of transcript variants using Entrez eutilities, you will have to use Entrez.elink to get a list of transcript accessions for the gene of interest and then use Entrez.efetch to fetch the sequence(s). Just to give you an idea, you can use Entrez Direct for this as follows:
elink -db gene -id 682 -target nuccore -name gene_nuccore_refseqrna \
  | efetch -format fasta

You can find more Entrez Link descriptions here. Note, there is a gene_fasta option for efetch but it still returns only the transcript sequence (i.e., no intron sequence).
If you need the full gene sequence (including the intronic regions), you can use the commands mentioned here.
Alternatively, you can use NCBI Datasets for this. You can use the web interface, command line tool or the python API. This page with the how-to articles has the information you seek.


In GFF3, annotating more than one protein-coding gene (i.e. polycistronic) contained in a eukaryotic mRNA

The GFF3 file format specification doesn't care where annotations are described in the file, although it can make some things easier if the annotations are ordered by the start position:

1   T1  gene    3631    4605    .   +   .   ID=ATNG01010
1   T1  mRNA    3631    4605    .   +   .   ID=ATNG01010.1;Parent=ATNG01010
1   T1  exon    3631    3913    .   +   .   ID=ATNG01010:exon:1;Parent=ATNG01010.1
1   T1  CDS     3750    3904    .   -   0   ID=ATRV01011:CDS:1;Parent=ATRV01011.0
1   T1  CDS     3860    3913    .   +   0   ID=ATNG01010:CDS:1;Parent=ATNG01010.1
1   T1  exon    3996    4276    .   +   .   ID=ATNG01010:exon:2;Parent=ATNG01010.1
1   T1  CDS     3996    4260    .   +   2   ID=ATNG01010:CDS:2;Parent=ATNG01010.1
1   T1  exon    4486    4605    .   +   .   ID=ATNG01010:exon:3;Parent=ATNG01010.1

There are many examples in well-studied genomes of overlapping genes, so the GFF3 format needs to be flexible enough to handle these. Different genes can be encoded across the same region, but on different strands. But different genes can also be encoded on the same strand, sometimes sharing an exon.


Assuming that it is the postprocessing script calculate_methylation_frequency.py that you are using, user Wouter de Coster's comment is correct that the methylation event is filtered out at the frequency estimation stage if it has a low LLR. It will have  a low LLR if it is considered not confident enough to report. This means not only LLR > T * N, but also the statistical threshold LLR > 0 (see below).
Specifically, here is the code in the script:

def update_call_stats(key, num_called_cpg_sites, is_methylated, sequence):
    if key not in sites:
        sites[key] = SiteStats(num_called_cpg_sites, sequence)

    sites[key].num_reads += 1
    sites[key].called_sites += num_called_cpg_sites
    if is_methylated > 0:
        sites[key].called_sites_methylated += num_called_cpg_sites

That is the function that decides whether to add a site to output, which depends on is_methylated, which is set here:
    # [skip several lines of code]
        num_sites = int(record['num_motifs'])
        llr = float(record['log_lik_ratio'])

        # Skip ambiguous call
        if abs(llr) < args.call_threshold * num_sites:
            # [MP NOTE: SKIP if LLR < T * N]
            continue
        sequence = record['sequence']

        # [MP NOTE: SKIP IF LLR <= 0]
        is_methylated = llr > 0

        # [skip several lines of code]

        update_call_stats(key, num_sites, is_methylated, sequence)

It's a little messy to read because the sites dictionary is a global outside of update_call_stats().
So overall, it seems like if you want sites included that are below this threshold, you will have to either:

use an earlier stage of processing (output of nanopolish call-methylation I believe according to the docs) or
alter the script to change the heuristic thresholds.


If you have issue with memory and dealing with large object, maybe data.table is the way to go (https://github.com/Rdatatable/data.table/wiki):

Let's take two dummy data.frame:

dfA <- data.frame(Chrom = rep("chr1",150),
                  Pos = 1:150,
                  ref = rep("A",150))
dfB <- data.frame(Chrom = rep("chr1",10),
                  Pos = 11:20,
                  Vaf = rep("VAF", 10))

The big data.frame dfA:

> head(dfA)
  Chrom Pos ref
1  chr1   1   A
2  chr1   2   A
3  chr1   3   A
4  chr1   4   A
5  chr1   5   A
6  chr1   6   A

and the small dataset:

> dfB
   Chrom Pos Vaf
1   chr1  11 VAF
2   chr1  12 VAF
3   chr1  13 VAF
4   chr1  14 VAF
5   chr1  15 VAF
6   chr1  16 VAF
7   chr1  17 VAF
8   chr1  18 VAF
9   chr1  19 VAF
10  chr1  20 VAF

Now, using data.table that is more efficient and less memory consuming for the manipulation of dataframe, you can do:

library(data.table)
setDT(dfA)
setDT(dfB)
dfA[dfB, on = c("Chrom","Pos"), vaf := i.Vaf]

And now, your data.frame dfA looks like:

> dfA[5:25]
    Chrom Pos ref  vaf
 1:  chr1   5   A <NA>
 2:  chr1   6   A <NA>
 3:  chr1   7   A <NA>
 4:  chr1   8   A <NA>
 5:  chr1   9   A <NA>
 6:  chr1  10   A <NA>
 7:  chr1  11   A  VAF
 8:  chr1  12   A  VAF
 9:  chr1  13   A  VAF
10:  chr1  14   A  VAF
11:  chr1  15   A  VAF
12:  chr1  16   A  VAF
13:  chr1  17   A  VAF
14:  chr1  18   A  VAF
15:  chr1  19   A  VAF
16:  chr1  20   A  VAF
17:  chr1  21   A <NA>
18:  chr1  22   A <NA>
19:  chr1  23   A <NA>
20:  chr1  24   A <NA>
21:  chr1  25   A <NA>
    Chrom Pos ref  vaf

As you see, we incorporate the column Vaf from the small dataframe into the big dataframe.

With your example, I will first rename the colnames of your small dataframe in order to have the perfect match with the big dataframe:

colnames(data) <- gsub("Chrom", "Chromosome", colnames(data))

then,

setDT(data1)
setDT(data)
data1[data, on = c("Chromosome","Position"), i_TumorVAF_WU := i.VAF]

Does it look what you are trying to obtain ?


not gnuparallel but nextflow. I don't know macs and what is 'first.bam' and 'second.bam' is not clear not me but I hope you'll get the idea.

The nextflow:

Channel.fromPath(params.bams).splitCsv(header: false,sep:'\t',strip:true).map{T->file(T[0])}.into{bams1;bams2}

boolean isControl(file) {
    return file.name.endsWith("_Input_sorted.bam_rem.bam") || file.name.endsWith("_IgG_sorted.bam_rem.bam ");
    }

bams1.filter{T->isControl(T)}.set{controls}
bams2.filter{T->!isControl(T)}.set{experiment}

process mac {
    maxForks 14
    tag "${bam1} vs ${bam2}"
    input:
        set bam1 , bam2  from experiment.combine(controls)
    output:
        file("output.peaks") into mac_peak
    script:
    """
    # check syntax please
    macs2 callpeak -t {bam1.toRealPath()}  -c ${bam2.toRealPath()} -g hs -n out-model -q 0.01 --outdir output.peaks
    """
    }

invoke:

find /dir1/dir2 -type f -name "*.bam" > input.list

nextflow run --bams input.list workflow.nf


Calculating celltype proportion changes between conditions in single-cell data

This publication used a negative binomial regression model of the counts of cells at each time-point.

From the methods section: "For each cell type, we model the number of lineage-labelled cells detected in each analysed mouse as a random count variable using a negative binomial distribution. The frequency of detection is modelled by using the natural log of the total number of cells of that type profiled in a given mouse as an offset. The time point of each mouse (0, 30 or 60 days post tamoxifen) is provided as a covariate. The negative binomial model was fit using the R command ‘glm.nb’ from the ‘MASS’ package. The P value for the significance of the change in labelled fraction size between time-points was assessed using a likelihood-ratio test, computing using the R function ‘anova’."

The method is available here as an R script.

The input x is a dataframe with the following columns:

x$variable # celltype
x$value # cell count
x$batch # condition+donor
x$timepoint # condition (e.g. adult, fetal)

Montoro, D.T., Haber, A.L., Biton, M. et al. A revised airway epithelial hierarchy includes CFTR-expressing ionocytes. Nature 560, 319–324 (2018) doi:10.1038/s41586-018-0393-7


Questions regarding Nanopore sequencing analysis

1. The quickest way to know if it's a single or multi fast5 is the file name. If it mentions channels and read numbers (e.g. ..._read_9_ch_471_..., then it's a single read file. Alternatively, the multi-read files have only read IDs at the top level (i.e. read_<32-character-string-with-additional-dashes>). Example:

$ h5ls <read_file_name>.fast5 | head
read_0115573a-042a-42a3-b4e2-f0338b4a6c66 Group
read_015b870d-8a03-4f22-a73e-fb9f833b94c7 Group
read_0167393a-dede-4a27-91af-494622b63261 Group
read_01e658de-853b-4f58-93ad-87fbd735254c Group
read_0202acc2-33f4-4d29-9dac-c7f80c758c86 Group
read_022d68f0-6f8e-4124-bfb0-11262efa4a67 Group
read_0258a3c9-4ffe-42ab-a42d-89466fcefbf5 Group
read_032da6db-e84a-43ec-8c25-bcaead992066 Group
read_03383754-be43-420b-87fc-60c2eb833963 Group
read_036b4824-bb0a-42dd-ade0-29d89112ff07 Group

2. ONT provide a single_to_multi program for converting to the multi-read format. Details can be found on this github page. Here's a quick example of how it could be run:

$ single_to_multi_fast5 --input_path fast5 --save_path fast5_multi \
    --filename_base multi_file_myRunID --batch_size 4000 --recursive

3a. For multi-fast5 files, the kit/cell attributes are in the context_tags group of the file (because a multi-fast5 file could contain fast5 elements from multiple different sequencing runs with different run parameters):

$ h5dump -g /read_<readID>/context_tags multi_file_myRunID_0.fast5 | \
  grep -A 11 -e 'sequencing_kit' -e 'flowcell_type'
   ATTRIBUTE "flowcell_type" {
      DATATYPE  H5T_STRING {
         STRSIZE 11;
         STRPAD H5T_STR_NULLTERM;
         CSET H5T_CSET_ASCII;
         CTYPE H5T_C_S1;
      }
      DATASPACE  SCALAR
      DATA {
      (0): "flo-min106"
      }
   }
--
   ATTRIBUTE "sequencing_kit" {
      DATATYPE  H5T_STRING {
         STRSIZE 11;
         STRPAD H5T_STR_NULLTERM;
         CSET H5T_CSET_ASCII;
         CTYPE H5T_C_S1;
      }
      DATASPACE  SCALAR
      DATA {
      (0): "sqk-pbk004"
      }
   }

3b. For single fast5 files (at least, in the most recent version of them), they are found in the UniqueGlobalKey group:

$ h5dump -g /UniqueGlobalKey single_file_read_X_ch_Y_strand.fast5 | \
  grep -A 11 -e 'sequencing_kit' -e 'flowcell_type'
      ATTRIBUTE "flowcell_type" {
         DATATYPE  H5T_STRING {
            STRSIZE 11;
            STRPAD H5T_STR_NULLTERM;
            CSET H5T_CSET_ASCII;
            CTYPE H5T_C_S1;
         }
         DATASPACE  SCALAR
         DATA {
         (0): "flo-min106"
         }
      }
--
      ATTRIBUTE "sequencing_kit" {
         DATATYPE  H5T_STRING {
            STRSIZE 11;
            STRPAD H5T_STR_NULLTERM;
            CSET H5T_CSET_ASCII;
            CTYPE H5T_C_S1;
         }
         DATASPACE  SCALAR
         DATA {
         (0): "sqk-pcs108"
         }
      }

[note that this file structure may change in the future]


Illustrate a 3D visualisation of the three main PCs using plot3d() package in R?

Regarding your question on how to decide which PCs to choose for plotting, prcomp is ordering PCs by their proportions in the variance. For example, using the dataset iris:

iris_pc <- prcomp(iris[,1:4], scale  = TRUE, center = TRUE)

and you get

> summary(iris_pc)
Importance of components:
                          PC1    PC2     PC3     PC4
Standard deviation     1.7084 0.9560 0.38309 0.14393
Proportion of Variance 0.7296 0.2285 0.03669 0.00518
Cumulative Proportion  0.7296 0.9581 0.99482 1.00000

BTW, using the princomp instead of prcomp gives similar results:

iris_pc2 <- princomp(iris[,1:4], cor=TRUE, scores=TRUE)

> summary(iris_pc)
Importance of components:
                          Comp.1    Comp.2     Comp.3      Comp.4
Standard deviation     1.7083611 0.9560494 0.38308860 0.143926497
Proportion of Variance 0.7296245 0.2285076 0.03668922 0.005178709
Cumulative Proportion  0.7296245 0.9581321 0.99482129 1.000000000

So, this means if you want to represent the main 3 PCs, you have to pick the first three PCs.

However, based on your dataset, if you take a look at the output of prcomp, you will see:

pc = prcomp(t(pcaData), scale=T, center=T)

> summary(pc)
Importance of components:
                          PC1    PC2     PC3     PC4     PC5     PC6     PC7    PC8     PC9    PC10   PC11    PC12
Standard deviation     9.3653 7.9806 5.46230 5.02379 3.95681 3.36495 3.21719 3.1042 2.98449 2.73437 2.6529 2.56841
Proportion of Variance 0.1993 0.1447 0.06781 0.05736 0.03558 0.02573 0.02352 0.0219 0.02024 0.01699 0.0160 0.01499

You can see that the first three PCs count for only 41% of the total variance. This can be explain by the diversity of your samples (coming from different cancers for example). So it means that your graphic representation won't be perfect to separate clusters.

Here, the graphic representation of the cumulative proportion of variances accross PC:

plot(pc, type = "l")

Next for the plotting and color attribution, I get a lot of trouble to extract the values from your Subgroups row.... for some reasons, I was not able to convert them into a character vector... weird ... anyway, I re-defined groups based on the code you put in patients ids:

patients = as.character(colnames(data))
grp = c(rep("C1",length(grep(".C1_",patients))),
        rep("C2",length(grep(".C2_",patients))),
        rep("C3",length(grep(".C3_",patients))),
        rep("C4",length(grep(".C4_",patients))),
        rep("C5",length(grep(".C5_",patients))),
        rep("C6",length(grep(".C6_",patients))))

For plotting results of PCA, you can use plot3d from rgl package for example. It gives you an interactive plot.

library(rgl)
colors = c("red","yellow", "green","cyan", "blue","purple")
plot3d(pc$x[,1:3], col =colors[as.factor(grp)], main = "pc / prcomp")
legend3d("topright",legend = unique(grp),col = colors, pch = 16, cex = 2)

And here the snapshot of your plot:

Hope that it helps you to figure it out the solution of your problem


parsimony and maximum likelihood tree comparison in R

Can i even compare these two methods?

See if this paper helps:
https://www.nature.com/articles/nature02917


The following sed command will do the trick.

sed -e '/^[^>]/ s/[a-z]/\-/g' in.fna > out.fna

Lines starting with > will be escaped and the lowercase characters will be replaced with - .


Seurat VlnPlot presenting expression of multiple genes in a single cluster

To do so one workaround it to have your data in "long format" and then use the column that holds the "gene names" as the x variable while plotting.

You can use FetchData() to extract data from a Seurat object. VlnPlot's default is the data slot (of the active assay if using Seurat v3 I suppose). And you can specify which cells and genes to retrieve.

selected_cells <- names(panc8$celltype[panc8$celltype == "gamma"])
data <- FetchData(panc8,
                    vars = c('FZD9','CTNNB1','APC'),
                    cells = selected_cells ,
                    slot = "data")
> head(data)
           FZD9   CTNNB1      APC
D101_5        0 0.000000 0.000000
D101_43       0 2.007853 1.001958
D101_93       0 0.000000 0.000000

Melt() would transform your data into "long format". By not specifying any arguments, we push for all of the info in the three variables to be gathered in two columns:

long_data <- melt(data)
No id variables; using all as measure variables
> head(long_data)
  variable value
1     FZD9     0
2     FZD9     0
3     FZD9     0

> tail(long_data)
     variable value
1870      APC     0
1871      APC     0
1872      APC     0

ggplot2 is used to add "violin" and "jitter" layers. You can customize the output to look (exactly) like the VlnPlot() output.

ggplot(long_data,
       aes(x = variable, y = value)) +
  geom_violin() +
  geom_jitter(size = 0.1)

And here is the same graph generated by VlnPlot():

There are no "violins" as the counts are almost entirely zeros. And see how different the ranges are in the y axes of the VlnPlot. So the "tweak" I have presented here would only work for genes that are expressed at similar levels / similar ranges.


The difference between the two depends on to whom you talk ;)

You are right: both refer to one base difference from the sequence. Usually the term "mutation" is used if the change has an impact on the phenotype. The "P" in "SNP" means "polymorphism". A lot of people use this term if more than 1% of the individual in a given population have this variant.

My recommendation: only use "mutation" if you are sure that your variant has an important impact on the phenotype. Only use "SNP" if you are sure that it is a commonly occurring variant. In any other case just use "variant". It is slowly becoming more common to talk about "SNV" (Single Nucleotide Variant) than about "SNP".


What is the current state-of-the-art in assembling hybrid transcriptomes?

My recommendation would be to run Trinity using the --long_reads option, which allows you to provide an error-corrected fastq file for anchoring reads to a transcript isoform. The reads are then clustered per-isoform, and assembled, in a similar fashion to what's done with genome-guided Trinity:
https://github.com/trinityrnaseq/trinityrnaseq/wiki/Running-Trinity
Trinity --seqType fq --max_memory 50G --long_reads corrected_pacbio_reads.fa --left reads_1.fq  --right reads_2.fq --CPU 6

Disclaimer: I have contributed code to Trinity at some time in the distant past


Approximate matches with bedtools intersect

You can use bedtools intersect but first use bcftools norm.  bcftools norm will left-align and normalize indels to get both vcf files consistent with each other.  This should resolve the matching if they are truly matching variants.


Error given while trying to index a BAM file with Samtools Index - NO COOR?

By now this will be of little interest to the original poster, but for future reference:
$ samtools sort -n -T ~/[etc] -o sample${r}bowtiemapping_sorted.bam [etc]
$ samtools index sample${r}bowtiemapping_sorted.bam
[E::hts_idx_push] NO_COOR reads not in a single block at the end 1652 -1.

This error message indicates (admittedly somewhat obscurely) that the reads are not sorted by coordinate — in particular, there are mapped reads appearing after the unplaced unmapped reads (those with RNAME *) that should be at the very end of the file.
This is because samtools sort -n has been used to sort the reads by name instead. Remove -n to sort by position, which is what is needed to prepare a BAM file for indexing with samtools index.


Why the t-test for a specific gene shows different value compared to differential analysis?

There is no reason your t-test should reproduce edgeR. In fact, edgeR exists because t-test is inappropriate.

edgeR does the tests by pooling information from all genes, because with the low number of replicates your t-test doesn't give sufficient statistical power.

What you need to do is: check visually your gene and make a decision yourself based on graphical results and edgeR results if your gene shows biologically significant.


Database for germline copy number variations in cancer

You can find CNV (and other structural variant) information in:

DECIPHER

  DECIPHER is used by the clinical community to share and compare phenotypic and genotypic data. The DECIPHER database contains data from 29603 patients who have given consent for broad data-sharing; DECIPHER also supports more limited sharing via consortia. Have a look at the numbers.

  Anyone can browse publicly-available patient data on DECIPHER and request to be put in contact with the responsible clinician. Why? Because sharing benefits everyone.

ExAC

  The Exome Aggregation Consortium (ExAC) is a coalition of investigators seeking to aggregate and harmonize exome sequencing data from a wide variety of large-scale sequencing projects, and to make summary data available for the wider scientific community.

  The data set provided on this website spans 60,706 unrelated individuals sequenced as part of various disease-specific and population genetic studies. The ExAC Principal Investigators and groups that have contributed data to the current release are listed here.

ExAC has included CNV calls since August 2016.
ClinVar

  ClinVar is a freely accessible, public archive of reports of the relationships among human variations and phenotypes, with supporting evidence. ClinVar thus facilitates access to and communication about the relationships asserted between human variation and observed health status, and the history of that interpretation. ClinVar processes submissions reporting variants found in patient samples, assertions made regarding their clinical significance, information about the submitter, and other supporting data.

Database of Genomic Variants

  The objective of the Database of Genomic Variants is to provide a comprehensive summary of structural variation in the human genome. We define structural variation as genomic alterations that involve segments of DNA that are larger than 50bp. The content of the database is only representing structural variation identified in healthy control samples.

  The Database of Genomic Variants provides a useful catalog of control data for studies aiming to correlate genomic variation with phenotypic data. The database is continuously updated with new data from peer reviewed research studies. We always welcome suggestions and comments regarding the database from the research community.

Information from all of these is also displayed in the Structural Variant browser of VarSome*. While you can't search by patient status, so you can't use VarSome to know the clinical history of any sample with the CNV, you can at least look up a specific gene and see the structural variants listed for it. For example, for TP53, you would see:

* Disclaimer: I work for the company behind VarSome.


What is the difference between samtools, bamtools, picard, sambamba and biobambam?

The obvious answer is that different people wrote them. It's fairly common in bioinformatics for people with a computer science background to get frustrated with existing tools and create their own alternative tool (rather than improving an existing tool). Over time, tools with similar initial aims will have popular functionality implemented in them (and eventually have bugs fixed), such that it matters less which particular tool is used for common methods.

Here's my impression of the tools:

samtools -- originally written by Heng Li (who also wrote BWA). The people who now work on samtools also maintain the alignment file format specification for SAM, BAM, and CRAM, so any new file format features are likely to be implemented in samtools first.
bamtools -- this looks like it was written by Derek Barnett, Erik Garrison, Gabor Marth, Michael Stromberg to mirror the samtools toolkit, but using C++ instead of C
picard -- Java tools written by the Broad Institute for manipulating BAM/SAM files. Being written in Java makes it easier to port to other operating systems, so it may work better on Windows systems. I'm more familiar with picard being used at a filtering level (e.g. removing PCR duplicates), and for statistical analysis, but it links in with the Java HTS library from samtools, so probably shares a lot of the functionality.
sambamba -- a GPL2-licensed toolkit written in the D programming language (presumably by Artem Tarasov and Pjotr Prins). I haven't used it (and don't know people who have used it), but the github page suggests "For almost 5 years the main advantage over samtools was parallelized BAM reading. Finally in March 2017 samtools 1.4 was released, reaching parity on this."
biobambam -- written by German Tischler in C++. I also have no experience with this toolkit. This seems to have some multithreading capability, but is otherwise similar to other toolkits.


How GFF3 attributes (9th column) varies from one gene prediction algorithm to another

The first place to start is the GFF3 specification. This is the official word on what is and is not allowed in a GFF3 file. For example, users can define arbitrary attribute keys, so long as they do not begin with an uppercase letter (these are reserved for "official" use).

But your question doesn't seem to be about what is allowed, but what is commonly used. I have a question for you: which gene predictors are YOU using? Or are you using gene annotations produced by others?

Here are a few examples.

NCBI RefSeq

If you're using GFF3 from NCBI, chances are it looks something like this.

NC_007070.3     Gnomon  gene    329235  331223  .       +       .       ID=gene14;Dbxref=BEEBASE:GB42168,GeneID:551678;Name=LOC551678;gbkey=Gene;gene=LOC551678;gene_biotype=protein_coding
NC_007070.3     Gnomon  mRNA    329235  331223  .       +       .       ID=rna27;Parent=gene14;Dbxref=GeneID:551678,Genbank:XM_624067.4,BEEBASE:GB42168;Name=XM_624067.4;gbkey=mRNA;gene=LOC551678;model_evidence=Supporting evidence includes similarity to: 66 ESTs%2C 24 Proteins%2C and 99%25 coverage of the annotated genomic feature by RNAseq alignments%2C including 113 samples with support for all annotated introns;product=receptor expression-enhancing protein 5-like%2C transcript variant X1;transcript_id=XM_624067.4
NC_007070.3     Gnomon  exon    329235  329459  .       +       .       ID=id117;Parent=rna27;Dbxref=GeneID:551678,Genbank:XM_624067.4,BEEBASE:GB42168;gbkey=mRNA;gene=LOC551678;product=receptor expression-enhancing protein 5-like%2C transcript variant X1;transcript_id=XM_624067.4
NC_007070.3     Gnomon  exon    329850  330082  .       +       .       ID=id118;Parent=rna27;Dbxref=GeneID:551678,Genbank:XM_624067.4,BEEBASE:GB42168;gbkey=mRNA;gene=LOC551678;product=receptor expression-enhancing protein 5-like%2C transcript variant X1;transcript_id=XM_624067.4
NC_007070.3     Gnomon  exon    330166  330301  .       +       .       ID=id119;Parent=rna27;Dbxref=GeneID:551678,Genbank:XM_624067.4,BEEBASE:GB42168;gbkey=mRNA;gene=LOC551678;product=receptor expression-enhancing protein 5-like%2C transcript variant X1;transcript_id=XM_624067.4
NC_007070.3     Gnomon  exon    330376  331223  .       +       .       ID=id120;Parent=rna27;Dbxref=GeneID:551678,Genbank:XM_624067.4,BEEBASE:GB42168;gbkey=mRNA;gene=LOC551678;product=receptor expression-enhancing protein 5-like%2C transcript variant X1;transcript_id=XM_624067.4
NC_007070.3     Gnomon  CDS     329333  329459  .       +       0       ID=cds8;Parent=rna27;Dbxref=GeneID:551678,Genbank:XP_624070.1,BEEBASE:GB42168;Name=XP_624070.1;gbkey=CDS;gene=LOC551678;product=receptor expression-enhancing protein 5-like isoform X1;protein_id=XP_624070.1
NC_007070.3     Gnomon  CDS     329850  330082  .       +       2       ID=cds8;Parent=rna27;Dbxref=GeneID:551678,Genbank:XP_624070.1,BEEBASE:GB42168;Name=XP_624070.1;gbkey=CDS;gene=LOC551678;product=receptor expression-enhancing protein 5-like isoform X1;protein_id=XP_624070.1
NC_007070.3     Gnomon  CDS     330166  330301  .       +       0       ID=cds8;Parent=rna27;Dbxref=GeneID:551678,Genbank:XP_624070.1,BEEBASE:GB42168;Name=XP_624070.1;gbkey=CDS;gene=LOC551678;product=receptor expression-enhancing protein 5-like isoform X1;protein_id=XP_624070.1
NC_007070.3     Gnomon  CDS     330376  330416  .       +       2       ID=cds8;Parent=rna27;Dbxref=GeneID:551678,Genbank:XP_624070.1,BEEBASE:GB42168;Name=XP_624070.1;gbkey=CDS;gene=LOC551678;product=receptor expression-enhancing protein 5-like isoform X1;protein_id=XP_624070.1

MAKER

The MAKER annotation workflow (paper, software) is a pretty commonly used gene annotation tool, and produces GFF3 output like this.

scaffold_12     maker   gene    652527  655343  .       +       .       ID=maker-scaffold_12-augustus-gene-0.959;Name=maker-scaffold_12-augustus-gene-0.959
scaffold_12     maker   mRNA    652527  655343  .       +       .       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1;Parent=maker-scaffold_12-augustus-gene-0.959;Name=maker-scaffold_12-augustus-gene-0.959-mRNA-1;_AED=0.24;_eAED=0.18;_QI=0|0|0|0.66|0.5|0.33|3|0|218
scaffold_12     maker   exon    652527  652817  .       +       .       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:exon:1203;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   exon    654877  655170  .       +       .       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:exon:1204;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   exon    655272  655343  .       +       .       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:exon:1205;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   CDS     652527  652817  .       +       0       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:cds;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   CDS     654877  655170  .       +       0       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:cds;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   CDS     655272  655343  .       +       0       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:cds;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   gene    941547  943897  .       +       .       ID=snap-scaffold_12-processed-gene-0.851;Name=snap-scaffold_12-processed-gene-0.851
scaffold_12     maker   mRNA    941547  943897  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1;Parent=snap-scaffold_12-processed-gene-0.851;Name=snap-scaffold_12-processed-gene-0.851-mRNA-1;_AED=0.22;_eAED=0.22;_QI=0|0|0.25|0.25|1|1|4|661|95
scaffold_12     maker   exon    941547  941631  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:exon:1206;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   exon    942343  942367  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:exon:1207;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   exon    942780  942920  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:exon:1208;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   exon    943200  943897  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:exon:1209;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   CDS     941547  941631  .       +       0       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:cds;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   CDS     942343  942367  .       +       2       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:cds;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   CDS     942780  942920  .       +       1       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:cds;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   CDS     943200  943236  .       +       1       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:cds;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   three_prime_UTR 943237  943897  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:three_prime_utr;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1

GeneMark, SNAP, Augustus

If we look at the output of several common ab initio gene prediction tools, none of them actually produces GFF3 by default. Here is some output from GeneMark (paper, software).

543   1   +  Initial         799051     799097         47          1 2 - -
543   2   +  Terminal        799236     799266         31          3 3 - -

544   3   -  Terminal        802357     802514        158          3 2 - -
544   2   -  Internal        802607     802685         79          1 1 - -
544   1   -  Initial         802829     802843         15          3 1 - -

Here is some output from SNAP (paper, software).

Einit   8230    8239    +       9.329   0       1       0       scaffold_12-snap.4
Exon    8848    8869    +       2.484   2       2       2       scaffold_12-snap.4
Exon    10121   10208   +       15.302  1       0       2       scaffold_12-snap.4
Exon    11361   11420   +       5.969   0       0       2       scaffold_12-snap.4
Exon    11471   11535   +       0.921   0       2       1       scaffold_12-snap.4
Eterm   12169   12187   +       18.163  1       0       1       scaffold_12-snap.4
Einit   14569   14668   +       -6.918  0       1       0       scaffold_12-snap.5
Exon    15029   15203   +       -1.023  2       2       0       scaffold_12-snap.5
Exon    16171   16307   +       -9.230  1       1       1       scaffold_12-snap.5
Eterm   16667   16698   +       9.829   2       0       0       scaffold_12-snap.5
Einit   17809   17898   -       -5.390  0       0       0       scaffold_12-snap.6
Exon    17299   17350   -       6.978   0       1       2       scaffold_12-snap.6
Eterm   17029   17129   -       -2.857  2       0       0       scaffold_12-snap.6

Both of these formats appear to be tab-delimited, but that's where the similarity with GFF3 ends.

The output of Augustus (paper, software) looks like GTF (a variant of the GFF format) once you lose all the lines beginning with a # symbol, but Augustus also has a GFF3 output mode.

# Predicted genes for sequence number 1 on both strands
# start gene g1
scaffold_12     AUGUSTUS        gene    2841    3400    1       +       .       g1
scaffold_12     AUGUSTUS        transcript      2841    3400    1       +       .       g1.t1
scaffold_12     AUGUSTUS        start_codon     2841    2843    .       +       0       transcript_id "g1.t1"; gene_id "g1";
scaffold_12     AUGUSTUS        intron  3027    3097    1       +       .       transcript_id "g1.t1"; gene_id "g1";
scaffold_12     AUGUSTUS        CDS     2841    3026    1       +       0       transcript_id "g1.t1"; gene_id "g1";
scaffold_12     AUGUSTUS        CDS     3098    3400    1       +       0       transcript_id "g1.t1"; gene_id "g1";
scaffold_12     AUGUSTUS        stop_codon      3398    3400    .       +       0       transcript_id "g1.t1"; gene_id "g1";
# protein sequence = [MAIKNAEHDLRVIVDAIEGLGLKVAPHKTEAMAFPASALCGRRGAAPPKIRLGGSSILVGSRSRWYISGHSENSSKSP
# RTEGKETTPLQQRDPLDAPLWVSGVVAHCCGGPEGQEGCPGLAAQGSDQGVLRIRDGLLCGYDGCGDHRPRPSDSSAGGGLCRP]
# end gene g1
###
# start gene g2
scaffold_12     AUGUSTUS        gene    4712    15229   0.21    +       .       g2
scaffold_12     AUGUSTUS        transcript      4712    15229   0.21    +       .       g2.t1
scaffold_12     AUGUSTUS        start_codon     4712    4714    .       +       0       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        intron  4858    5591    0.74    +       .       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        intron  5686    15028   0.23    +       .       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        CDS     4712    4857    0.72    +       0       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        CDS     5592    5685    0.42    +       1       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        CDS     15029   15229   0.43    +       0       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        stop_codon      15227   15229   .       +       0       transcript_id "g2.t1"; gene_id "g2";
# protein sequence = [MGRNSHRSCCVVNCKITSAKSDCKFYKFPTAKWKINQRKMWVAAVKRQKYIKDEISHAETQTEITEVTGATKVNYANK
# KYICLLFVRTYVLRMLVDVALSNLRFSLFGIRKSLEIFGQSEKADQTRWRLPSCEMEWIESRKGKMRE]
# end gene g2
###

Parsing attributes

So...how should YOU handle attributes in GFF3's 9th column? That depends a lot on what you want to do with the data.

The most important attributes are ID and Parent, which are used to define relationships between features and subfeatures. (These relationships implicitly define a directed acyclic graph of features, although most GFF3 parsers don't directly support traversal of this graph.)

But not only do you have to handle the pre-defined attributes discussed in the GFF3 specification, you also have to be able to handle any number of arbitrary attributes whose keys you may not know beforehand.

The simplest way to handle this would be to parse the attribute column into a dictionary of key/value pairs. Once it's in this form, it's trivial to see what attributes are there and how to access them.

For example, if we have a feature with the following attributes...

ID=mRNA42;Parent=gene19;integrity=0.95;foo=bar

...we would want it in a dictionary like so.

attributes = {
    'ID': 'mrna42',
    'Parent': 'gene19',
    'integrity': '0.95',
    'foo': 'bar',
}

Python code to parse that might look something like this.

attributes = dict()
for keyvaluepair in attributestring.split(';'):
    for key, value in keyvaluepair.split('='):
        attributes[key] = value

One thing to consider is that it's possible for an attribute to have multiple values (separated by commas). It's not commonly used, but it is valid and something that should be handled or at least checked for.


processing fastq files using cellranger on linux

this has been resolved. I downloaded the windows file and this caused trouble but when downloading the file within the Linux, it seems to be working.


What's the difference between VCF spec versions 4.1 and 4.2?

This is easy to check, you can download both specs in .tex format and do diff.

Changes to the v4.2 compared to v4.1:

Information field format: adding source and version as recommended fields.
INFO field can have one value for each possible allele (code R).
For all of the ##INFO, ##FORMAT, ##FILTER, and ##ALT metainformation, extra fields can be included after the default fields.
Alternate base (ALT) can include *: missing due to a upstream deletion.
Quality scores, a sentence removed: High QUAL scores indicate high confidence calls. Although traditionally people use integer phred scores, this field is permitted to be a floating point to enable higher resolution for low confidence calls if desired.
Examples changed a bit.


What tools can I use for a bacterial core/pan genome pipeline?

LS-BSR should be able to give you what you're looking for. See the article.

After you run the primary analysis, there is a simple documented workflow for splitting the pangenome into core and accessory, based on a user defined threshold. I'm the developer, so can help if you run into any problems.


What process and input data is required for a cellranger reference transcriptome?

Your problem is caused by using the transcriptome fasta file rather than the genome fasta file. You've already given it transcriptome information with genome_genes.filtered.gtf, it needs the genomic sequence to go along with that.

As an aside, your GTF filtering command is unlikely to be useful unless you have a Gencode or Ensembl GTF for human. Have a look at the gene_biotypes in your file and decide for yourself what (if anything) needs to be filtered. For example, I made a cellranger reference last night for dm3 (the old fruit fly genome) with:

cellranger mkgtf annotation.gtf genes.gtf \
    --attribute=gene_biotype:protein_coding \
    --attribute=gene_biotype:ncRNA \
    --attribute=gene_biotype:snRNA \
    --attribute=gene_biotype:snoRNA \
    --attribute=gene_biotype:pre_miRNA

That should be customized for relevance. It's quite possible that your GTF doesn't even have gene_biotype attributes, that's specific to Ensembl I think.


First off,

Don’t use RPKMs.

They are truly deprecated because they’re confusing once it comes to paired-end reads. If anything, use FPKMs, which are mathematically the same but use a more correct name (do we count paired reads separately? No, we count fragments).

Even better, use TPM (= transcripts per million), or an appropriate cross-library normalisation method. TMP is defined as:

$$
\text{TPM}_\color{orchid}i =
    {\color{dodgerblue}{\frac{x_\color{orchid}i}{{l_\text{eff}}_\color{orchid}i}}}
    \cdot
    \frac{1}{\sum_\color{tomato}j \color{dodgerblue}{\frac{x_\color{tomato}j}{{l_\text{eff}}_\color{tomato}j}}}
    \cdot
    \color{darkcyan}{10^6}
$$

where

$\color{orchid}i$: transcript index,
$x_i$: transcript raw count,
$\color{tomato}j$ iterates over all (known) transcripts,
$\color{dodgerblue}{\frac{x_k}{{l_\text{eff}}_k}}$: rate of fragment coverage per nucleobase ($l_\text{eff}$ being the effective length),
$\color{darkcyan}{10^6}$: scaling factor (= “per millions”).

That said, FPKM an be calculated in R as follows. Note that most of the calculation happens in log transformed number space, to avoid numerical instability:

fpkm = function (counts, effective_lengths) {
    exp(log(counts) - log(effective_lengths) - log(sum(counts)) + log(1E9))
}

Here, the effective length is the transcript length minus the mean fragment length plus 1; that is, all the possible positions of an average fragment inside the transcript, which equals the number of all distinct fragments that can be sampled from a transcript.

This function handles one library at a time. I (and others) argue that this is the way functions should be written. If you want to apply the code to multiple libraries, nothing is easier using ‹dplyr›:

tidy_expression = tidy_expression %>%
    group_by(Sample) %>%
    mutate(FPKM = fpkm(Count, col_data$Lengths))

However, the data in the question isn’t in tidy data format, so we first need to transform it accordingly using ‹tidyr›:

tidy_expression = expression %>%
    gather(Sample, Count)

This equation fails if all your counts are zero; instead of zeros you will get a vector of NaNs. You might want to account for that.

And I mentioned that TPMs are superior, so here’s their function as well:

tpm = function (counts, effective_lengths) {
    rate = log(counts) - log(effective_lengths)
    exp(rate - log(sum(exp(rate))) + log(1E6))
}


Is the number of contigs and scaffolds an artifact of genome assembly and rather arbitrary, or is it canonical for each species?

This is purely arbitrary and an artifact of how well the assembly turned out. Better quality assemblies with have fewer and longer contigs.


What is cellranger doing in comparison to other methods?

cellranger mkfastq is not necessary anymore.  It used to be that the cellranger software wanted the reads to be interleaved, and you could use cellranger to do that for you if you couldn't do it yourself.  Newer versions of cellranger will take the fastq files just like Illumina's bcl2fastq makes them.

Cellranger count aligns the reads, filters away duplicates based on UMIs, tries to figure out which cell barcodes really captured cells, and tells you for each cell barcode how many reads hit each gene.  You can go through the .bams, and see which reads were counted for what gene and what barcode, based on the flags and other tags in the .bam

If you want another pipeline to compare to, you can use the McCarroll lab's DropSeq pipeline, which is a bunch of little scripts, around the output of STAR aligner (just like 10XGenomics uses).  I have put the same data sets through both 10X and Dropseq, and they return about the same thing.


Measure the purifying selection for certain taxa along a phylogeny

The following is not the clearest explanation of dN/dS humankind has seen, but I hope you get the gist and core idea of likelihood is absolutely correct. I am wondering if you are thinking of a different test to the way I would do this calculation. This ain't a text book clarity answer.
Anyway what you are proposing is a likelihood ratio test (LRT) for  selection based on the parameterisation of omega. Your idea isn't quite correct ... its not dN/dS value you are looking at to assess significance, its the likelihood. This point is critical and is the most important difference in our thinking on this test.
LRT works because its a log likelihood and subtracting two log likelihood is the same as dividing them (I forget the theory - but you subtract the likelihood to work out the p-value via Chi squared). There's some theory about why the Chi-squared distribution works but I can't remember it.
You need to work out the degrees of freedom which is the parameter differences between 'free branch model' and the neutral model and this value is LOW. If the neutral model is the "global model", then its 2 degrees of freedom - IN YOUR TREE. This is because the SP you are looking at form a single clade (monophyly). From memory you stick a #1 on the branch leading into the clade, ie. defining the node and thats your 'model' on the treefile. Here you have two dN/dS one for the whole tree excluding your clade and one for your taxa. This gives you a likelihood.
Note you could stick a #1 on each branch within the clade, i.e. the treefile. That means something different - like every branch was subject to ... in this case purifying selection ... not just the common ancestor.
You now work out the likelihood of the 'global' model... i.e. remove the #1 and make sure it gives you only 1 dN/dS value.
Subtract the likelihoods and multiple by two - if its LRT (super important).
You now need to work out the degrees of freedom. If its the same model but one has an additional node with a different dN/dS ... the df = 2. If you've dumped a #1 on every branch in the clade ... I still think df = 2 - you need to check.
If the difference in likelihood results in p < 0.05 you are correct the difference is significant ... whatever that difference might be.
If the neutral model is some different calculation of omega within the population genetics process (see here) then it will be the parameter differences between your model and the 'neutral' model.
To reiterate
I like simulations so much I type in reiterations ... So again the LRT usual will usually subtract one likelihood from another and THEM multiply the probability by two ... and compare the result on the Chi squared table for the difference in degrees of freedom. And again what you are looking at is not the value of dN/dS but the difference in LIKELIHOOD. That is very important. You need to be careful with the degrees of freedom on an LRT because they can introduce an additional subtraction constant, test dependings (its ages since I've done one so I can't remember).
If the test is significant and one dN/dS is x for the clade and y for the rest of the tree, then that differences is significant and if y is >1 (bit unlikely) and x < 1 then sure thats purifying selection. Usually clades evolve via positive selection rather than purifying selection.
From memory that is my understanding ... but if you have a neutral model of omega and a purifying model of omega ... yeah sort of, but usually thats for a 'global' model, i.e. no differences within the tree. In your case your interested in a specific clade. The idea is still the same its measured through likelihood.
Its so long since I've done this, but I do caution not to confuse two different tests.

I do understand Yang et al did an MCMC for dN/dS ... what I'm thinking of here is maximum likelihood. You'd need to watch for that carefully and thats what you are referring to when you say Chi-squared, because that is LRT 100%. LRT is maximum likelihood.


"megahit: command not found" after correct installation with miniconda

type conda info --envs.  Look at your list of environments.  Say it's base that has megahit in it, do conda activate base, then just do megahit.


Differential gene expression bias due to effect of an individual sample

Leave-one-out cross validation is an iterative method that leaves out one sample until each sample has been left out once. You can then identify the intersection of the differentially expressed genes from the sample space of repeated analyses to produce a set of genes that are consistently differentially expressed independent of individual sample biases.
more info here:
http://efavdb.com/leave-one-out-cross-validation/


Sequence alignment of extended IUPAC sequences

You can try using Biostrings, depends on how large is your reference, if it is small, you can just pairwiseAlignment with the nucleotideSubstitutionMatrix which supposedly :

library(Biostrings)
data(phiX174Phage)

x = DNAString("TCCAWRSCGCAGAAGT")

lapply(phiX174Phage,function(y)pairwiseAlignment(x,y,
substitutionMatrix = nucleotideSubstitutionMatrix()))

$Genbank
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92

$RF70s
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92

$SS78
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92

$Bull
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92

$G97
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92

$NEB03
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92

If it is much larger, you can convert the sequence of interest to a pwm and find it in the reference fasta, like shown in this example


How to isolate genes from whole genomes for phylogenetic tree analysis?

There are a lot of ways to do this. I suggest using Prokka/Roary to produce a core genome alignment. There's a useful tutorial on the Roary website:

for file in *
do
    prokka --kingdom Bacteria --outdir "${file%%.*}"  --genus Listeria --locustag "${file%%.*}" "$file"
    mv "${file%%.*}"/PROKKA_07052017.gff GFF/"${file%%.*}".gff # use current date
done
roary -f Alignment -e -n -v GFF/*.gff

Alignment/core_gene_alignment.aln can be used as input for phylogenetic analyses


When to account for the blacklisted genomic regions in ChIP-seq data analyses?

Aside: Cross-correlation is largely meaningless, regardless of what some of the ENCODE folks might argue. When we process our DEEP samples we don't even look at that value.

Regardless, if you're using SPP/phantomPeakQual for cross-correlation then note that it already removes the highest peaks from your dataset before computing the cross-correlation (in fact, it can remove most of the actual peaks too, which makes one further wonder what it's actually telling you). I don't know that this is actually documented anywhere, it's something I noticed when going through the code while pondering whether to implement it in deepTools. But at least it's ignoring these regions already :)

In general, it's most convenient to just remove peaks overlapping blacklisted regions. In an ideal world you'd filter out the blacklisted reads before peak calling, but (1) this is really inconvenient (more time and disk required) and (2) I've never seen an appreciable gain in peak calling performance. In theory at least you should be losing sensitivity right around blacklisted regions if you don't remove reads overlapping blacklisted regions, but you have to ask yourself whether you want to trust such peaks anyway. For other QC steps, at least with deepTools we provide a parameter with every tool to specify a BED file of blacklisted regions to skip.

As an aside, there are many fewer blacklisted regions in more recent genome builds (GRCh38 and GRCm38 at least), so this is less of an issue in general with them.


Transposable elements (TEs) are the classic example of rapidly evolving regions in eukaryotic genomes. These are often used as a proxy for neutrally evolving sequences, which are by definition not conserved/constrained. The UCSC genome browser includes various repeat tracks (e.g. RepeatMasker) that includes TEs.

Within protein coding regions, 4-fold degenerate sites are often used as a proxy for neutrally evolving sequences. These sites are where any nucleotide at this position leads to the same amino acid and therefore protein (synonymous mutations).

Potentially even more variable sequences would be those under positive selection, but there is not a single class of element that falls into this category i think.


Import a tab-separated file with differing numbers of elements in each row; prokka output

You don't really need to get around anything, your example can be loaded correctly with read.delim(). The rows with missing values are filled in with blanks. You may, however, prefer the readr package, which handles this a bit more elegantly (it'll tell you which values were missing and not do the annoying string->factor conversion).


Batch alignment of inconsistently named Fastq files

A very simple approach, but one which will work with the names you show, is to remove any occurrence of 00N (where N is any number) from right before the .gz. To do so, edit your script to:

for i in $(ls *.fastq*.gz | sed 's/00[0-9]\.gz/.gz/' | rev | cut -c 16- | rev | uniq); do
  bowtie2 -q -p20 -x /run/media/punit/data1/BowtieIndex/hg38 \
    -1 ${i}R1_001.fastq.gz \
    -2 ${i}R2_001.fastq.gz |
        samtools view -bS - >${i%}.bam"
done


The salmon output log contains a line with [info] Mapping rate which allows to calculate how many reads went mapped or unmapped. Is that sufficient? It should be in the logs folder of the directory salmon creates per sample.


How to use analysis output in my pipeline

I assume you are referring to the pepstats program from EMBOSS. The program has a -outfile flag that allows you to specify your desired name for the output file when you're not running the program interactively.


How to edit the headers of multiple fasta files from multiple folders? (recursively)

You can edit your for loop so as to include the paths to your files of interest.

The command below will give you the relative paths (relative to where you run the command, I think it is better to run this within your main directory) of your files of interest:

ls mpg*/mgm*.3/*.fna

You can use the output of this command (your files of interest) to iterate over via your for loop. Your code will look like:

for i in {...}; do
    sed -i.bak "s/^>/>${i}_/" i
done


How to do `bedtools intersection` using pandas alone?

Your setup:

import pandas as pd

dict1 = {0:['chr1','chr1','chr1','chr1','chr2'],
    1:[1, 100, 150, 900, 1], 2:[100, 200, 500, 950, 100],
    3:['feature1', 'feature2', 'feature3', 'feature4', 'feature4'],
    4:[0, 0, 0, 0, 0], 5:['+','+','-','+','+']}

df1 = pd.DataFrame(dict1)

print(df1)

##       0    1    2         3  4  5
## 0  chr1    1  100  feature1  0  +
## 1  chr1  100  200  feature2  0  +
## 2  chr1  150  500  feature3  0  -
## 3  chr1  900  950  feature4  0  +
## 4  chr2    1  100  feature4  0  +

dict2 = {0:['chr1','chr1'], 1:[155, 800], 2:[200, 901],
    3:['feature5', 'feature6'], 4:[0, 0], 5:['-','+']}

df2 = pd.DataFrame(dict2)

Now use pyranges (https://github.com/biocore-ntnu/pyranges):

# pip install pyranges
# conda install -c bioconda pyranges

import pyranges as pr
cols = "Chromosome Start End Name Score Strand".split()
# Add column names to the dfs so PyRanges knows which
# represent Chromosomes, Starts, Ends and optionally Strand

df1.columns = cols
df2.columns = cols

# create PyRanges-objects from the dfs
gr1, gr2 = pr.PyRanges(df1), pr.PyRanges(df2)

# intersect the two
gr = gr1.intersect(gr2)

# get the result as a   df
df = gr.df

Results:

print(gr)
# +--------------+-----------+-----------+------------+-----------+--------------+
# | Chromosome   |     Start |       End | Name       |     Score | Strand       |
# | (category)   |   (int32) |   (int32) | (object)   |   (int64) | (category)   |
# |--------------+-----------+-----------+------------+-----------+--------------|
# | chr1         |       155 |       200 | feature2   |         0 | +            |
# | chr1         |       900 |       901 | feature4   |         0 | +            |
# | chr1         |       155 |       200 | feature3   |         0 | -            |
# +--------------+-----------+-----------+------------+-----------+--------------+
# Stranded PyRanges object has 3 rows and 6 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome and Strand.
print(df)
#   Chromosome  Start  End      Name  Score Strand
# 0       chr1    155  200  feature2      0      +
# 1       chr1    900  901  feature4      0      +
# 2       chr1    155  200  feature3      0      -

See the docs for more info: https://biocore-ntnu.github.io/pyranges/


Classification (supervised learning) of expression data on pathway level

Have a look at the GSVA package. It allows to convert a matrix with genes x Samples to a pathways x Samples using several methods ssgsea, gage, gsva...

Afterwards you can use that matrix as input for differential expression of pathways or classification algorithms or whatever.

However it depends on the input of the "pathways" you give it. Make sure that the pathways are correct. Also consider that the pathway topology might change with conditions.


Why does my SplitDotPlotGG plot show differential expression of a gene but Seurat::FindMarkers() function does not return it?

You're showing groups 1 and 2 in the plot but comparing groups 0 and 1 in the differential expression test, so you need to instead set ident.1 = 1, ident.2 = 2 in FindMarkers.

The FindMarkers function by default only returns genes with >0.25 log fold-change expression, with >3 cells in each group. It's possible that these criteria are not being met for gene X. You can try changing these defaults, setting logfc.threshold = 0 to return all genes.


Remove Redundant Sequences from FASTA file in Python

Your code (as currently formatted) should generate the correct output, but it is not efficient. Make seen a set to improve these O(n) in calls to O(1) time complexity- this will speed it up a lot
seen = set()
records = []

for record in SeqIO.parse("b4r2.fasta", "fasta"):
    if record.seq not in seen:
        seen.add(record.seq)
        records.append(record)


Another 2 tree compare tools: ggtree and blatic 3


Generate table for total number of SV events per sample

So you have data in "long format" and want to convert it to "wide format". You can use spread() from tidyr for that:

require(tidyr)
# Assuming your input is called `foo`
spread(foo, Sample, value)


How to use the hmmsearch on prototypic repetitive sequences Repbase Update database?

RepBase explain what these files are on their website. It is usually a good idea to check the documentation when downloading files so you know what you're downloading.

Based on the list of files described on the linked page, there is no prirep.ref file, so I assume you mean plnrep.ref which is "Other Plants". In any case, again according to that page, the files with plant and/or fungal data (which is what I assume you want; you never said) are:

athrep.ref  Arabidopsis thaliana
oryrep.ref  Oryza sativa
dcotrep.ref Other Dicot Plants
mcotrep.ref Other Monocot Plants
grasrep.ref Other Grass Plants
plnrep.ref  Other Plants
fngrep.ref  Fungi

There's also chlrep.ref for the green alga Chlamydomonas reinhardtii, and diarep.ref for the diatom Thalassiosira pseudonana. I don't know whether you would want to consider either here.

In any case, since these are fasta files, you cna combine them by simply cating them together:

cat athrep.ref oryrep.ref dcotrep.ref mcotrep.ref grasrep.ref \
    plnrep.ref fngrep.ref > plantsAndFungi.ref

The, run your search on the plantsAndFungi.ref file.


How to get the longest fasta sequence including all possible switching isoforms of a gene out of isoforms

From your question it isn't clear if you are trying to identify the isoforms that have a very high or very low differential expression value, or if you are trying to correctly find an orthlog for your top hits using BLAST.

If your question is about isoform-level quantification, identifying DE at the isoform level is challenging and DESeq2 is only designed for gene-level analyses. That is why it is only outputting genes rather than a specific isoform.


How do I select a proper design matrix to use in DRIMSeq?

From your question I understand that you have two treatments in two mouse strains.
If you are interested in general treatment effects independet of the animal strain you would use the both variables in the design formula (similar to controlling for a batch effect).

meta.data <- data.frame(treatment = c("a", "a", "b", "b"), strain = c("x", "y", "x", "y"))
design_combined <- model.matrix(~ treatment + strain, data = meta.data)

If you want to make the comparisons for each strain separately, the easiest way will be to split your data set by strain and do a set up just for treatment. Alternatively, you can also merge the two columns into a new one in case you also want to also compare between the strains.

meta.data$combined <- paste0(meta.data$treatment, "_", meta.data$strain)
# refactor to change the reference level
meta.data$combined <- factor(meta.data$combined, levels = "bx", "by", "ax", "ay")
design_combined <- model.matrix(~ combined, data = meta.data)

To get the fold changes, you can extract the proportions for each transcript using

DRIMSeq::proportions(dataSet)

And from the proportions the log2 fold changes can be calculated with ease.

However, when I look at your plot something seems very off with your analysis right now. In all the examples I have seen and also in my own data the diamonds are representing the mean proportion of a transcript in a condition. See this example from the DRIMSeq manual on Bioconductor:

In your case there are several different means per group and if you look at the line to the very right the two treatments seem to be to be right on top of each other.
So I think that treatment and strain are somehow mixed in your setup and this is causing this kind of merged plot. If you use a design as outlined above you should get much better results. Additionally, you could maybe check that the meta data is correct.

Also in case you have not read it yet, the manual is a great resource:
https://bioconductor.org/packages/release/bioc/vignettes/DRIMSeq/inst/doc/DRIMSeq.pdf


I think you could use the MESH.db package in Bioconductor to annotate your genes. I have never used, but you should be able to do something like:

library("MeSH.db")
select(MeSH.db, keys = myGenes, keytype = "ENTREZ", columns = "MESH")

It depends on the MeSHDbi package from which you could retrieve other data form that database if needed


There are actually 2048 usable sequencing wells, hexagonally packed with four wells connected to the same sequencing sensor/channel via a multiplex (mux) selector. The combination of the mux and the channel number determines the physical location of the well. Unfortunately, the association between channel number and physical location is not obvious, and within each channel the muxes are also not in an obvious order, i.e. [3,4,1,2,2,1,4,3] for two adjacent channels.

Here's an "official" text description of the layout as it was for R7; it hasn't changed for the most recent flow cells (R9.5):

Channels 1-64 occur at the top of the chip (the other low channel numbers are at the bottom)
Channels order down the chip: 1-64, 449-512, 385-448, 321-384, 257-320, 193-256, 129-192, 65-128
Muxes run from left to right in the order: 3, 4, 1, 2, 2, 1, 4, 3

To confirm this layout, I pulled the channel/position lookup table definitions out of MinKNOW, and it seems to be the same as this:

Here's the code to produce this image (I haven't adjusted it to have a hexagonal structure):

#!/usr/bin/Rscript
library(rjson);

layout.fname <- "/opt/ONT/MinKNOW/Client/resources/app/resources/channels.json";

data.list <- fromJSON(file=layout.fname);

png("MinION_FC_Layout.png", width=3072, height=1536, pointsize=24);
par(mar=c(4,4,0.5,0.5), lwd=3);
data.mat <- matrix(unlist(data.list), ncol=2, byrow=TRUE);
plot(data.mat, pch=21,
     bg=(colorRampPalette(c("white","grey"))(2048))[1:2048],
     col=c("red","yellow","green","blue"),
     xlab="locX", ylab="locY", cex=2,
     xlim=c(7,76), ylim=c(0,33));
text(x=data.mat, labels=rep(1:512, each=4), cex=0.4);
legend("top",legend=c(1,2,3,4), fill=c("red","yellow","green","blue"),
       inset=0.05, title="MUX");
text(x=6, y=data.mat[data.mat[,1] == 8,2],
     labels = ceiling(which(data.mat[,1] == 8) / 4));
text(x=77, y=data.mat[data.mat[,1] == 75,2],
     labels = ceiling(which(data.mat[,1] == 75) / 4));
arrows(x0=c(8,75), x1=c(39,44), y0=33, length=0.2);
invisible(dev.off());


Subsetting of object existing of two samples

I think your issue is simply not assigning the newly-subset object in to a new object.

Assuming rownames(subset_DR1) is indeed a list of cell names

seuratObj_subset_dr1 <- SubsetData(seuratObj, cells=rownames(subset_DR1))

will work.


What is the index fastq file (sample_I*.fastq.gz) generated when demultiplexing Illumina paired-end runs?

tldr - The I*.fastq.gz file contains the read index sequences.
long explanation
Illumina uses a program called bcl2fastq to demultiplex sequencing runs.
This software takes a list of samples and their associated indices and uses those sequences to make one or more fastq files per sample, binned by one or two index sequences on either end of the sequencing molecule (i5 and i7 indices, see page 6 for HiSeq).
Illumina sequencing isn't perfect though, and there are sometimes errors in reading the index sequence. For example this index, CAGCCCAC can easily have read errors in the A sandwiched between many Cs: CAGCCCAC -> CAGCCCCC.
Instead of throwing out all of the reads with indexing sequencing errors, the bcl2fastq program includes reads that are discernibly derived from a true index for that sample as long as there is no overlap with another sample.
In the example provided above, the four indices for mysample were: ACATTACT, CAGCCCAC, GGCAATGG, TTTGGGTA. If we look at all of the index sequences in the Mysample_I1_001.fastq.gz file, we will see that sequences with the four correct indices are indeed the most prevalent, however there are also reads that have indices with sequencing errors derived from the correct indices.
bioawk -cfastx '{print($seq)}' Mysample_I1_001.fastq.gz | \
    sort | uniq -c | sort -k1 -nr | head

count      index       source
41362311   CAGCCCAC    True index
37209190   GGCAATGG    True index
36863213   ACATTACT    True index
33674467   TTTGGGTA    True index
 1140358   NAGCCCAC    CAGCCCAC
 1026099   NGCAATGG    GGCAATGG
 1016754   NCATTACT    ACATTACT
  933342   NTTGGGTA    TTTGGGTA
  119626   TTTGGGGA    TTTGGGTA
   98657   GTTGGGTA    TTTGGGTA
   96625   GGCAATGA    GGCAATGG


Usually the home of a program is where the executable is. In order to make the system aware you need to include it as a variable when the terminal is run.

So you need to add it to the file at ~/.bashrc of your computer (if you are using bash as shell [you can check it with echo $SHELL]). You need to add :

export RATT_HOME=path/where/RATT/is/installed/or/downloaded

As terdon says in a comment below:

  Note that echo $SHELL just prints the default login shell of the
  user, and not necessarily the shell that is currently running.  Also,
  ~/.bashrc is not read by macOS bash which starts login shells by
  default. In any case, this sort of global variable should be in
  ~/.profile, not .bashrc. For more of this sort of detail see
  Sequence of scripts sourced upon login and How to correctly add
  a path to PATH.


This is a bug affecting the curl compile, it's been reported here and also here, I suggest trying to compile from the newer release of the libcurl source rather than the .deb you linked.


Counting letters in phylip alignment columns with Biopython

using something like for col in column(alignment):

I may get the number of columns which have all the same letters (amino acids)

len(set(col)) == 1

The number of columns which have at least 2 different amino acids (with no dashes):

len(set(col)) > 1 and '-' not in col

The number of columns which have any single dash:

'-' in col

or if you really want columns with a single dash:

    col.count('-') == 1


Designing a lab NGS file database schema

The Global Alliance for Genomics and Health has been working on the issue of representing sequencing data and metadata for storage and sharing for quite some time, though with mixed results.  They do offer a model and API for storing NGS data in their GitHub repository, but it can be a bit of a pain to get a high-level view.  I am not sure if any better representation of this exists elsewhere.

I can say from personal experience (having built over a dozen genomic databases), there is no ideal data model and storage best practices.  Genomic data comes in many shapes and sizes, and your needs are going to vary from every other organization, so what works for one bioinformatics group won't necessarily work for you.  The best thing to do is design and implement a model that will cover all of the data types in your workflow and downstream analyses you might do with the data and metadata.


You could set up a public read-only bucket in an S3 account and put your files there. URLs would look like:

https://${username}.s3.amazonaws.com/public/${filename}

For samtools, the request for a BAM file's index should be built off that URL. For IGV, I think you can specify separate, custom URLs for the data and the index resources.

Public buckets are available to all, without need for credentials. So that might be an issue to think about, with respect to billing for web requests and data transfer.

I'm still not 100% certain that samtools/tabix support https requests reliably (I have a prior question on here about that), but that might be one approach.


What is the difference between a Bioinformatics pipeline and workflow?

Great question! Note that from a prescriptive standpoint, the terms pipeline and workflow don't have any strict or precise definitions. But it's still useful to take a descriptive standpoint and discuss how the terms are commonly used in the bioinformatics community.

But before talking about pipelines and workflows, it's helpful to talk about programs and scripts. A program or script typically implements a single data analysis task (or set of related tasks). Some examples include the following.

FastQC, a program that checks NGS reads for common quality issues
Trimmomatic, a program for cleaning NGS reads
salmon, a program for estimating transcript abundance from NGS reads
a custom R script that uses DESeq2 to perform differential expression analysis

A pipeline or a workflow refers to a particular kind of program or script that is intended primarily to combine other independent programs or scripts. For example, I might want to write an RNA-seq workflow that executes Trimmomatic, FastQC, salmon, and the R script using a single command. This is particularly useful if I have to run the same command many times, or if the commands take a long time to run. It's very inconvenient when you have to babysit your computer and wait for step 3 to finish so that you can launch step 4!

So when does a program become a pipeline? Honestly, there are no strict rules. In some cases it's clear: the 10-line Python script I wrote to split Fasta files is definitely NOT a pipeline, but the 200-line Python script I wrote that does nothing but invoke 6 other bioinformatics programs definitely IS a pipeline. There are a lot of tools that fall in the middle: they may require running multiple steps in a certain order, or implement their own processing but also delegate processing to other tools. Usually nobody worries too much about whether it's "correct" to call a particular tool a pipeline.

Finally, a workflow engine is the software used to actually execute your pipeline/workflow. As mentioned above, general-purpose scripting languages like Bash, Python, or Perl can be used to implement workflows. But there are other languages that are designed specifically for managing workflows. Perhaps the earliest and most popular of these is GNU Make, which was originally intended to help engineers coordinate software compilation but can be used for just about any workflow. More recently there has been a proliferation of tools intended to replace GNU Make for numerous languages in a variety of contexts. The most popular in bioinformatics seems to be Snakemake, which provides a nice balance of simplicity (through shell commands), flexibility (through configuration), and power-user support (through Python scripting). Build scripts written for these tools (i.e., a Makefile or Snakefile) are often called pipelines or workflows, and the workflow engine is the software that executes the workflow.

The workflow engines you listed above (such as argo) can certainly be used to coordinate bioinformatics workflows. Honestly though, these are aimed more at the broader tech industry: they involve not just workflow execution but also hardware and infrastructure coordination, and would require a level of engineering expertise/support not commonly available in a bioinformatics setting. This could change, however, as bioinformatics becomes more of a "big data" endeavor.

As a final note, I'll mention a few more relevant technologies that I wasn't able to fit above.

Docker: managing a consistent software environment across multiple (potentially dozens or hundreds) of computers; Singularity is Docker's less popular step-sister
Common Workflow Language (CWL): a generic language for declaring how each step of a workflow is executed, what inputs it needs, what outputs it creates, and approximately what resources (RAM, storage, CPU threads, etc.) are required to run it; designed to write workflows that can be run on a variety of workflow engines
Dockstore: a registry of bioinformatics workflows (heavy emphasis on genomics) that includes a Docker container and a CWL specification for each workflow
toil: a production-grade workflow engine used primarily for bioinformatics workflows


Differential Expression With 2 Treatments

If you're interested in looking at the differences between two treatments then you'll end up wanting to do both a direct contrast as well as the individual comparisons to baseline.

The direct contrast will give you the genes actually differentially expressed between the two conditions. In practice, you may want to filter this a bit so you only have genes differentially expressed vs. baseline in at least one condition (e.g., to get rid of genes only slightly higher due to T1 and slightly lower in T2, but not different enough in either case to be DE). Use a pretty lax p-value threshold (e.g., 0.1 or 0.2) for this.


Script to allow gene set enrichment analysis of 10x genomics data in R

Seruat will give you a list of genes which it thinks are upregulated in a particular cluster. Look at the functions that talk about marker genes - these functions basically do a DE analysis of the genes in one cluster compared to the others.

Then take that list and feed it to any standard GO analysis tool. Have a look at the topGO topKEGG and geneSetTest functions in the limma package, the GOStats package and the gsea. All should be suitable.  The GOSeq package is designed to compensate for gene length bias in RNA-seq experiments. As 10X only samples 3'-tags, there shouldn't be any gene-length bias in the data, so this shouldn't be an issue.


How to select only RNA with Hetero atoms from pdb file with python?

Solution I found: (c/p of body)

class RNASelect(Select):
    def accept_residue(self, residue):
        return 1 if ((residue.id[0] != 'H_ MG') and ((residue.get_resname() == '  A' or residue.get_resname() == '  U' or   residue.get_resname() == '  G' or  residue.get_resname() == '  C' or residue.get_resname() == '  T') or  (residue.id[0][0:2] == 'H_' and (residue.get_resname()[2] == 'U'  or residue.get_resname()[2] == 'A' or residue.get_resname()[2] == 'C' or residue.get_resname()[2] == 'G' or residue.get_resname()[2] == 'T'))))  else 0

It's not pretty but it do the job. The idea is to get everything that finish by a nucleotide. I figured out that hetatom's names of RNA chains always finish by a nucleotide letter so i just have to check everything that begins by 'H_' and ends by 'A,T,C,G,U'. Just remove 'H_ MG' and you're good.


Can a gff file be converted to a fasta file?

A gff is a file of annotation.  It generally doesn't include sequence information, so you can not inter-convert.  But the file on your page has a fasta entry at the bottom..  I'd just grab the bottom half of it, and if the sequence length looks right, that's likely what you want.


groupGO parameters explanation (ont and level)

As mentioned in the comments MF is molecular function, BP is biological process and CC is cellular compartment. These are the 3 domains of the ontologies. The level refers to the level in the ontology graph. In the example in that link, pigmentation would be level 1, pigmentation during development level 2, regulation of pigmentation during development is level 3 and so on. You can find more discussion of this in this post on biostars from a year ago.

Regarding the exact settings you should use, I don't personally find CC to be an informative domain. The first couple levels are typically too generic to be useful and they'll include too many genes anyway for you to likely pick up any changes. I'd start around level 3 and see how things go (in other programs you would get all of the levels and domains at once).


Why doesn't Biopython AlignIO.read() recognise the 'mauve' format?

You're using version 1.68 or older. Mauve support was added in 1.70.


How to filter info score post-imputation?

In simple GWAS setups, each SNP is analyzed independently.  In those cases you can filter out SNPs with poor INFO scores at any point.

For analyses that combine information across SNPs (for example FaSTLMM), I would recommend filtering out SNPs before running the step that combines information across SNPs.

So, if you filter right after running Impute2, you should be safe.


Which reference to use for read mapping for popular model organisms

I will answer one of the points – E.coli.

TL;DR Bacteria, and in particular E.coli, are highly variable and there is usually no single best assembly. Large scale WGS studies should come with multiple assemblies for individual monophyletic clusters.

Long answer: Whereas a single reference sequence can make sense for human or mice (to a certain extent), bacteria are not the case. The variability between them even within a single species can be enormous and various lineages can contain very different genes. This is the reason why biologists usually distinguish a core genome consisting of the genes shared among a vast majority (typically 95%) of individuals within a given phylogenetic clade, (i.e., of those genes which are essential for this clade – imagine a CPU and RAM in a computer), and an accessory genome consisting of the other, unnecessary, genes – imagine a joystick in our analogy). The size of the accessory genome can vary a lot for different bacteria. For instance, Escherichia coli is said to have an open accessory genome (i.e., much bigger than the core genome; imagine a Raspberry-PI and all the extension kits) whereas Chlamydia trachomatis has a closed accessory genome (i.e., small compared to the core genome; imagine an ATM).

For species having a small accessory genome, a single reference might be sufficient. However, if a core genome is large, we can use either pan-genomes, or multiple references. The word pan-genome tends to be used for different things: the computational biologists usually use it as a representation of all genomic content in a certain phylogenetic clade (which can be mathematically modeled using finite automata, hence as graphs) whereas some biologists define it as a concatenation of all (core and accessory) genes from the clade (without considering SNPs, etc.). Unfortunately, methods for graph references (especially for read mapping) are still not sufficiently developed, therefore, researchers still have to resort to reference sequence-based methods and use either concatenated genes or multiple reference sequences.

When particular bacterial species are studied using massive whole-genome sequencing of many isolates, researchers usually cluster the isolates to monophyletic sequence clusters and infer a special reference sequence for each of these clusters (see, e.g., this example for the pneumococcus). All of the inferred references should be then published together with the paper or in a separate data paper (see, e.g., the sequences from the previous example).

Disclaimer: I am not a biologist. If anything in this post is too simplified or incorrect, please, leave me a comment and I will try to fix it. I tried to make this text accessible for my younger self.


Number of edges in Graphical Analysis using Jupyter with Go Programming language

the dot file in the Git repo has 6646 edges, the pajek file in your link has 7182.
Chances are the latter has been updated since.


I cannot tell exactly from the post or the linked project what type of data these scripts are designed for, but I am going to assume that "Scripts involved in our workflow for detecting CNVs from WGS data using read depth-based methods" refers to short-read Illumina sequencing.

Simulations can be helpful, but they are usually imperfect. Nowadays we have long-read data available. I would suggest obtaining some long-read data for which short-read data is also available, doing CNV detection on both datasets independently, and comparing the results.

The problem with using long-read data is that any differences you find will be conflated with the biases of the technology. So, if you find a difference, then this approach will not tell you if the differences are due to the biases of the technologies, or due to bugs in the scripts. If you use simulations, then this conflation of biases remains.

Short of re-writing the pipeline, the only other suggestion I have is to borrow from software testing. Try constructing a few minimal examples that might also test some corner cases and see if the scripts behave exactly as expected.


How do I find identical sequences in a FASTA file?

The trick would be to swap the key in the dictionary to be the sequence itself.

Also I would recommend using a different separator that "_" since that is what the current ids have so that you can easily separate the individual ids from the concatenated id. I used a pipe "|" in this example. Also I just manually wrote the FASTA output instead of using Biopython to do it.

from Bio import SeqIO
from collections import defaultdict

dedup_records = defaultdict(list)
for record in SeqIO.parse("test.fasta", "fasta"):
    # Use the sequence as the key and then have a list of id's as the value
    dedup_records[str(record.seq)].append(record.id)
with open("Output.fasta", 'w') as output:
    for seq, ids in dedup_records.items():
        # Join the ids and write them out as the fasta
        output.write(">{}\n".format('|'.join(ids)))
        output.write(seq + "\n")


How to manipulate protein interaction network from String database in R?

I think the easiest way is to download the graph using STRINGdb.

library(STRINGdb)
string_db <- STRINGdb$new(version="10", species=9606,
                          score_threshold=400, input_directory="" )
full.graph <- string_db$get_graph()

Now you can use igraph, to manipulate the graph. Let's assume you want to take 200 proteins with the highest degree, i.e. number of edges they have.

library(igraph)

# see how many proteins do you have
vcount(full.graph)

# find top 200 proteins with the highest degree
top.degree.verticies <- names(tail(sort(degree(full.graph)), 200))

# extract the relevant subgraph
top.subgraph <- induced_subgraph(full.graph, top.degree.verticies)

# count the number of proteins in it
vcount(top.subgraph)

How to get disease specific genes?

There's no GO annotation for cancer or Alzheimer's disease. It is out of scope of the GO consortium.

What you can do, you can either take KEGG Pathways annotation, or manually select list of relevant GO-terms. Or acquire the list from one of the papers. For example annotation term 05200 corresponds to the cancer KEGG pathway. You can easily retrieve proteins associated with the annotation:

cancer.pathway.proteins <-
    string_db$get_term_proteins('05200')$STRING_id

And then perform subgraphing as described above.

Alternatively you can try to get an enrichment score for an every gene given it's neighbors (the way enrichment is shown on the string-db website). Then you can keep only those having top enrichment scores. Probably get_ppi_enrichment_full or get_ppi_enrichment functions will help you to do that.


Does it make sense to compute values such as averages, sums, and differences (e.g. the average fold change for treatment1)? And yes, wary that doing things on log fold-change is not a thing, but how about just fold change?

This is problematic, if you want to compare treatments then compute the actual fold change and its significance. You can't really base anything off of the fold-change (log or otherwise) alone.

does it also make sense to compute something like Euclidean distances for clustering fold change data

It's not unheard of to do this. The problem is you lose all information about how reliable those fold-change estimates are. This is something maintained when you cluster sample-level information (e.g., TPM or robust log transformed data). Not that using variance-stabilized data for clustering will work better, since it tamps down on how much highly-expressed genes influence your clustering (assuming you don't want that to occur).
Have a look at the DESeq2 vignette for some nice examples of what good practices are.


Downloading Fasta for each chromosome of a complete genome

It is common that reference genomes will be published only as a complete genome, without separate files for each chromosome. Fortunately, a google search for "split fasta by chromosome" will give you many options for downloading a complete genome Fasta file and splitting it by chromosome. The top answer on this biostar thread looks promising.

pip install pyfaidx
pyfaidx -x sequences.fa

Alternatively, I typically use a homebrew Python script. If there are unlocalized/unplaced scaffolds in the file, this will create a bunch of files with ugly unexpected names. These are easy to ignore or delete.

# usage: ./split.py < genome.fasta
from Bio import SeqIO
import sys
reader = SeqIO.parse(sys.stdin, 'fasta')
for record in reader:
    outfile = 'chr{chrid}.fasta'.format(chrid=record.id)
    with open(outfile, 'w') as outstream:
        count = SeqIO.write(record, outstream, 'fasta')

With a little bit of searching, you will find a wealth of helpful information on blogs, forums, and library documentation.


Why is there no GT-entry in my .vcf file from bacteria

The representation for a variant present in a haploid genome should be 1. A heterozygote in a diploid genome would be 0/1 (or 0|1 or 1|0  if phased). So, to add a 1 as the genotype, you would need to first insert a GT to the beginning of each INFO line and then the 1 at the beginning of the last field. But, this should happen only on lines that are not headers (do not start with #). This should do it:
$ awk -F'\t' -v OFS="\t" '{ if(/^#/){print}else{$9="GT:"$9;$10="1:"$10;print}}' file.vcf
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Variants:_LGC19-XL01_S63_L001_R_001_(trimmed)
CFC381_K12_Bw25113  360287  .   GGGAAT  AATTCC  1952.76 .   NS=1;VF=0.380;SB=0.533;SB50=0.60;SB65=1.0;TYPE=Substitution;AVQUAL=30   GT:DP:AO    1:231:89

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

Now, in this specific case, you could also do something like this:
$ sed -i.bak '/^[^#]/{s/\<DP:AO\t/GT:DP:AO\t1:/}' file.vcf

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


How can I edit a specific FASTQ read in place, given the read ID?

You'll first need to determine the appropriate line number, which you can do with grep -m1 -nw "@31027" foo.fastq. After that, note that you can provide a line number to sed:

sed -i '123456s/CAAAAGTCACTCA/CAAAAGTCACTTA/' foo.fastq

That will do the replacement only on line 123456 and edit the file in place (the -i option).


What is this 5-column sequencing format from 2009?

This is an early Solexa/Illumina sequencer format. The columns are the identifier to location on the flowcell.  I believe the "." was the original placeholder for an unread base, which has been replaced by and "N" in current Illumina sequencing output.
From on http://www.crg.eu/en/content/processing-and-analysis-illumina-sequencing-data.

seq.txt (Gerald)
4 1 23 1566 ACCGCTCTCGTGCTCGTCGCTGCGTTGAGGCTTGCG
The first 4 columns are an ID that provides location details of the cluster in the flow cell, followed by the read sequence. Seq.txt file size typically ranges between 250 Megabytes and 1 Gigabyte. One file is generated per lane.


How to concatename file names in Snakemake?

You can use join to convert the list of input file into a string.
I think this would work as expected.:

rule run_prog:
    input: "a", "b"
    output: "output"
    run: f"prog {','.join(input)}"


Add segregating sites to branches of a tree

The segregating sites in parsimony are formally called steps. Any number of programs will do this but I don't really use R, particularly for phylogeny, You can of course simply look at the raw treefile, all the data is there (described in the last paragraph). If look at the options for genetic distance, in parsimony this will be an integer and again the technical name is steps and is the segreating sites you want.
Python's ETE3 package will work (but you'd need to pipe it through Biopython), as well as specific standalone packages Figtree, the steps or 'genetic distance' on the branches on the branches option is in the check box. The heavy weight program for tracking step associated phylogeny using parsimony is Mesquite and its famous grandparent is called MacClade.
I suspect its a nomenclature thing, the options probably talk about genetic distance because few people do parsimony phylogeny now. Thats the option you want  and you will know because it will be >>1, whilst distance and model-based methods will be between 0 to 1. From a tree viewing it doesn't matter if the genetic distance is 0 to 1 or 1 and beyond because its all proportional. Its just an option thing you are missing. Thus if there is an option, e.g. to include bootstrap support on the tree, there'll be an option for genetic distance.
Just make sure you don't confuse the number of steps per branch with a consensus statistic, because there're completely different but might give similar looking outputs. Again its easily checked against the raw treefile because the "steps" will be immediately adajacent to the taxa name for the "twigs" (terminal branches) in the tree.


getSequence has only been enabled for the main Ensembl (vertebrates) biomaRt. This is all at the Bioconductor end. You'll need to use a getBM instead, for example:

os_genes <- getBM(attributes =c("ensembl_transcript_id", "peptide"), filters = "ensembl_transcript_id", values = osj_id, mart = osj)


How to get % similarity between strains and mutation files

left_lineagelist = [x.split('_')[-1].split('|')[0]
                          for x in left_lineagelist]
        right_lineagelist = set([x.split('_')[-1].split('|')[0]
                          for x in right_lineagelist])

Allows for extraction of 6 digit EPI, provided the file has had sequences removed prior; as such:
for line in lines:
    if line[0] == '>':
        print(line[1:])


Is it possible to identify cells that are expressing two or more genes in Seurat?

The counts stored in the Seurat object are: raw counts (seuratobject@raw.data), the log + normalized counts (seuratobject@data), and the scaled counts (seuratobject@scale.data). FeaturePlot() plots the log + normalized counts.

In order to identify double-positive cells, you need to identify cells that express a gene (i.e. positive for a gene) and that is not trivial: you need a cutoff value for each gene, and also to distinguish between true and dropout zero counts.

Provided you solve these, you can retrieve your cells by:

MYH11.cutoff <- 1
ACTA2.cutoff <- 1
length(which(FetchData(seuratobject, vars.all='MYH11') > MYH11.cutoff))
length(which(FetchData(seuratobject, vars.all='ACTA2') > ACTA2.cutoff))
length(which(FetchData(seuratobject, vars.all='MYH11') > MYH11.cutoff & FetchData(seuratobject, vars.all='ACTA2') > ACTA2.cutoff))


Problem with classification model of genomic data: every machine learing model predicts wrongly almost always the same subset of dataset

You could be seeing the effect of SNPs on haplotype-driven conditions (i.e. where a combination of multiple variants predicts the class, rather than a single variant).

Do you have any method to determine the age of particular variants (e.g. via an LD plot to find inconsistencies between variants and the surrounding block structure), so that you can see if that correlates with the likelihood of error?


You could take a look at kraken2 It should fit your purpose.


Smith-Waterman traceback with large gaps

This question is based on a misunderstanding. The Smith-Waterman algorithm calculates the best local segment. You have to start at the $max(H_{i,j})$ point and run the traceback as described. It will (most likely) return only a local segment that, however, has the maximum value in terms of alignment score.
To get a global alignment, another algorithm has to be used. Needleman-Wunsch would be an example.
The difference between local and global is that a local alignment only alignes a subsequence of the provided sequences whereas a global alignment gives the best alignment of the whole sequences.


DESeqDataSetFromTximport invalid rownames length

I got exactly the same problem not long ago (see my other question what I tried to do). We used the same function!

My code ran fine on my HPC but the same code/data failed on my laptop. I later found out I had to upgrade my R to the latest version (3.5.1) and upgrade all my bioconductor packages to the latest. My HPC was on a recent bioconductor docker, so there was no problem.

Please make sure everything in your environment latest (R, R packages, Bioconductor).