When working with Bioconductor and the curatedOvarianData dataset, what are the unit values of the expression data? I am wondering how to interpret values inside the gene expression levels contained within Bioconductor. For example, if we have the following commands via R, we get:

source("https://bioconductor.org/biocLite.R")
biocLite("curatedOvarianData")
library(curatedOvarianData)
data(TCGA_eset)
head(exprs(TCGA_eset)["CXCL12",])

TCGA.20.0987 TCGA.23.1031 TCGA.24.0979 TCGA.23.1117 TCGA.23.1021 TCGA.04.1337
    5.437923     6.307250     4.193748     9.071293     5.101532     6.697770

The above calls out gene expression levels for the gene CXCL12. However, I am wondering how to interpret the values 5.437923, 6.307250, etc.

I know they are expression levels, but do they come with units? Do larger values mean greater cancer possibility?

These values represent the expression of a particular gene in each sample. While expression of an immune system gene like CXLC12 might be associated with the aggressiveness of a cancer, that's not going to be the case for the majority of genes.

Looking at the numbers, these look like they are probably log$_2$-transformed values. In order to know for sure, plotting a density distribution of the dataset is a good idea:

plot(density(exprs(TCGA_eset)));

If the values are in the range of around 0..15, then it'll be log$_2$-transformed values, whereas if they're in the range of 0..100000, it'll be "raw" expression values. The log-transformed values are preferred for data analysis because gene expression tends to be close to a normal distribution in log-space.

These numbers are best treated as relative counts (i.e. the most appropriate unit is "expression"), as variation in the preparation of the sample (or the sequencer, or the microarray) can influence the results that come out the other end. Converting them to something more meaningful (e.g. "transcripts per cell") requires normalising against a known reference sample.


Looper: missing argument error While trying to run looper I am getting the following error:

  Error (missing attribute): 'sandro_rna_seq' requires sample attribute '/scratch/nv4e/sandro_looper/pipeline/jyvo_experiment-metadata.yaml' for option/argument '-E'
  Not submitted: Required attribute(s) missing for pipeline arguments string

In my pipeline_interface.yaml file I have the following:

arguments:
  "-L": "IU"
  "-T": "gencode_mouse_m13"
  "-E": "/scratch/nv4e/sandro_looper/pipeline/jyvo_experiment-
metadata.yaml"

The file "/scratch/nv4e/sandro_looper/pipeline/jyvo_experiment-metadata.yaml" is definitely present but it still gives an error...

Any suggestions would be greatly appreciated.

You are using the argument values incorrectly. Argument values are sample attributes, not full character paths. Hence, "IU" is a sample attribute, "gencode_mouse_m13" is a sample attribute, so why would you think "/scratch..." would be treated any differently?

The pipeline interface maps sample attributes to command-line arguments. It is modulated by sample (different values per sample). What you're trying to do is pass a constant string to an argument.

For this, there is a pipeline_args section of your project config. Please read the docs on project config.


Possibility to save output blastn table in memory using biopython Is there a possibility, using biopython to save output table of blastn in memory, not in file on hard drive, for process it using pandas for example and then delete it from memory?

I mean, is there a possibility to redirect output of "out" argument of this function into memory? In pandas DataFrame for example.

NcbiblastnCommandline(query="fasta1.fasta", subject="fasta2.fasta", out="output_table.csv", outfmt=6)

I always saved output tables in file before, but now I need to align one set of small fasta files against another set of small fasta files and I think that the best way is avoid writing of many small tables in hard drive because after recieving one table I have to read and process this table immediately and after that I don't need this table anymore.

You can specify send to stdout using out='-' in the Biopython wrapper.

from Bio.Blast.Applications import NcbiblastnCommandline
import pandas as pd

cline = NcbiblastnCommandline(query='seq.fna', subject='seq2.fna', outfmt=6, out='-')
output = cline()[0].strip()

rows = [line.split() for line in output.splitlines()]
cols = ['qseqid', 'sseqid', 'pident', 'length',
        'mismatch', 'gapopen', 'qstart', 'qend',
        'sstart', 'send', 'evalue', 'bitscore']

data_types = {'pident': float, 'length': int, 'mismatch': int,
              'gapopen': int, 'qstart': int, 'qend': int,
              'sstart': int, 'send': int, 'evalue': float,
              'bitscore': float}

df = pd.DataFrame(rows, columns=cols).astype(data_types)


Is there any way to align ChIP-seq reads to telomeres? I know that telomeres are highly repeated sequences, but is there any way to retain any reads that map to these regions (on HG38)?

I recently managed to find some protein binding to centromeres, which are also mainly repeats. Therefore I wondered if there was any way to pick up signal on the telomeres.

My team and I are investigating binding of a protein complex to DNA and the experimental team reported that they saw it a bit on telomeres, so we were wondering if there was any way to measure this through bioinformatics. If not then it's no problem, we would just like to know.

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


Tool to show DNA sequence and allowing upload of own graph data Background

We want to be able to load (or request) data for a genome including the sequence and gene annotation (bacteria). Then, we want to load our own annotation which should be displayed as a line plot:

position score
1        5
2        6
etc.

Where each position corresponds to the position in the genome (i.e. FASTA).
Something like this:

Question

Is there a tool that can load a genome sequence, gene annotation, and a position ~ score file, which can be displayed as a line plot?

You can use UCSC genome browser, which supports custom genome sequence and annotation. Specifically, one needs to create a Track Hub to hold genomes and annotations.

Here is an example:

  http://genome.ucsc.edu/cgi-bin/hgTracks?hubUrl=http://zlab-trackhub.umassmed.edu/yfu/insects/hub.txt&db=hub_129615_cl_v0.8.2

And here is the corresponding hub.txt file holding the metadata, and URLs to actual fasta, bigWig and bigBed files.


Modeling number of reads mapped to a gene I am looking for a probability distribution of a number of reads mapped to a particular gene in metagenomic sequencing (NGS, shotgun, likely illumina).
Naively one could model it via a binomial (or Poisson) distribution as:
$$
P_n = {N\choose n}\pi^n(1-\pi)^{N-n},
$$
where

$n$ is the number of reads mapped to the gene
$N$ is the sequencing depth (total number of reads)
$\pi\propto l$ is the probability of read mapping to the particular gene, proportioal to the gene length (possibly adjusted for the read length)

This is however complicated by such intricacies of sequencing as:

finite length of the fragments, which can be longer or shorter than the gene length
paired-end sequencing, which means one has to be careful when counting evens where only one or both reads map to a gene
something else?

Is there an established model for this case? Could you recommend reading on the subject?

This is quite a similar question to this one, but I notice that you're specifically asking about a probability distribution rather than differential expression. That sounds like gene dispersion and/or shot noise; a good reference for me on that has been the documentation for DESeq2:
http://www.bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
The gene dispersion issue is discussed in more detail in their 2014 paper:
https://doi.org/10.1186/s13059-014-0550-8

We model read counts K ij as following a negative binomial distribution (sometimes also called a gamma-Poisson distribution) with mean μ ij and dispersion α i . The mean is taken as a quantity q ij, proportional to the concentration of cDNA fragments from the gene in the sample, scaled by a normalization factor s ij , i.e., μ ij =s ij q ij . For many applications, the same constant s j can be used for all genes in a sample, which then accounts for differences in sequencing depth between samples. To estimate these size factors, the DESeq2 package offers the median-of-ratios method already used in DESeq. However, it can be advantageous to calculate gene-specific normalization factors s ij to account for further sources of technical biases such as differing dependence on GC content, gene length or the like, using published methods, and these can be supplied instead.


Installing HiC-Pro issues I am setting up a HiChIPseq pipeline, but installing software necessary for the HiC side of things has been unbecoming.
Both my PI and I have tried to install Hi-C Pro onto our Mac Pro that we use as our compute server (though in a couple of weeks we hope to migrate to an EPYC 64 Core computer server of our own!).
There are three approaches that I have taken to install Hi-C Pro.

From Source

To install from source, you need to supply the paths to its dependencies. Believing we did a good job, we run it and find problems accessing bx-python.
No matter how many times we update pip and try to install bx-python and specify that we want it installed for python2.

From singularity

This is the one I tried. From what I understand this creates a VM with the necessary dependencies installed, and allows you to use the program without building. However, it appears to abort when conda tries to install numpy.
+ conda install -y -c anaconda numpy=1.16.3
Collecting package metadata (current_repodata.json): done
Solving environment: failed with initial frozen solve. Retrying with flexible solve.
Collecting package metadata (repodata.json): / Killed
ABORT: Aborting with RETVAL=255
Cleaning up...

I have seen the failed solving environment bug many times, but have no idea how to solve it. I have tried updating conda, but that doesn't solve it.

Install Someone's Premade Conda Environment

Running conda.sh give me this:
Solving environment: failed

ResolvePackageNotFound:
<Lots of packages...>

Sorry if this is sloppy.

I could install HiC-Pro using conda when creating a separate environment (which I always do for side-projects), and then pull the conda build of HiC-Pro:
conda create --name HiCPro python=2.7
conda activate HiCPro
conda install -c davebx hicpro

Not sure what conda.sh is that you describe, sometimes tricky to get random code from GitHub running. Give it a go with a standard conda install as I described.


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

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

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

Thank you in advance.

This may not be a viable answer for all experimental contexts, but you can compare expression profiles of individual cells against datasets that contain known reference cell types (migratory and not). We have recently published a web-based tool called CIPR (cluster identity predictor) that helps with cluster annotations in scRNAseq experiments. There are various reference datasets available, but you can also provide a custom reference by preparing data from databases like GEO.
You can read the manuscript for more information about how the algorithm works, but in a nutshell, it compares differentially expressed genes or global expression patterns within your data to known reference datasets and creates quick visual outputs. If your reference contains both migratory and non-migratory cell types, maybe CIPR can help differentiate them. For other identity prediction tools, check out this article is a good resource that compares different cluster calling algorithms.


What are the right parameters to trim a small RNA transcriptome with trimmomatic? I'm having some problems with finding the right parameters to trim my small RNA Illumina reads (51 nt long) with Trimmomatic.

Before trimming, one of the samples (21M reads) looks like this:

So for my understanding, it's quite good and the last two issues can be easily solved with adapter trimming.

I then tried to trim it with the command

java -jar /opt/Trimmomatic-0.38/trimmomatic-0.38.jar SE -threads 24 -phred33 FemaleMito3.fastq FemaleMito3_trimmed.fastq ILLUMINACLIP:adapters.ultimate.fa:2:30:10 AVGQUAL:25 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:18

but I got the following result:

So now I have a problem with the sequence quality per tile (why?), an overall drop of the quality at the end of the reads and a warning with the sequence length distribution.

Then I tried trimming just the adapters with this command:

java -jar /opt/Trimmomatic-0.38/trimmomatic-0.38.jar SE -threads 24 -phred33 FemaleMito1.fastq FemaleMito1_noadapters.fastq ILLUMINACLIP:adapters.ultimate.fa:2:30:10 AVGQUAL:25

but again I get some new issues:

Again the per tile sequence quality issue, and now a warning about the GC content.

I didn't find anywhere some precise instruction on how to trim a small RNA
library with Trimmomatic. How should I properly trim my small RNA library with Trimmomatic?

The command where you trim with adapters and by quality is perfectly fine. That FastQC isn't perfectly happy is expected. It's a tool made for whole genome sequencing QC, you SHOULD see a number of "fail"s with any RNA-seq protocol. Steps that should always fail in RNA-seq:

per-base sequence content (there's "random priming" performed that isn't completely random, so the first ~8 bases show a typical bias pattern).
per-sequence GC content (particularly in sRNA-seq you're enriching for a subset of genes, they won't have some nice GC distribution and that's perfectly OK).
Sequence duplication level (if this doesn't fail then everything is uniformly expressed...which would be quite unlikely)
Overrepresented sequences (as above, though this will be worse for sRNA-seq since reads will often cover the entire transcript, so there's even less sequence diversity)


Filtering bases based on phred qualities with pysam Is there a way to filter bases in BAM files based on phred quallities through python's pysam ?

I have a code here that

Takes the nucleobases per position from a BAM file using pysam's pileup function
Saves it in ReverseList and ForwardList based on both strands (i.e forward and reverse), I want to reject those bases that have phred quality below 25 to be not stored in the ForwardList and ReverseList lists so that they are not used for further analysis.

 samfile = pysam.Samfile( filename, "rb" )
 ReverseList = [''] * lenref
 ForwardList = [''] * lenref
 for pileupcolumn in samfile.pileup() :
     for pileupread in pileupcolumn.pileups:
         if (pileupread.alignment.mapping_quality <= 15):
               continue
         if not pileupread.is_del and not pileupread.is_refskip:
               if pileupread.alignment.is_reverse: #negative
                      ReverseList[pileupcolumn.pos] += pileupread.alignment.query_sequence[pileupread.query_position]
               else:
                      ForwardList[pileupcolumn.pos] += pileupread.alignment.query_sequence[pileupread.query_position]

 samfile.close()

Where lenref = 16569(length of mitochondrial genome) and filename is the name of BAM file. I want to filter based of phred qualities of bases.

The PileupRead object as a query_position attribute, which you can use for this:

for pileupcolumn in samfile.pileup() :
    for pileupread in pileupcolumn.pileups:
         if (pileupread.alignment.mapping_quality <= 15):
             continue
         if not pileupread.is_del and not pileupread.is_refskip:
             if pileupread.alignment.query_qualities[pileupread.query_position] < 10:
                 # Skip entries with base phred scores < 10
                 continue
             if pileupread.alignment.is_reverse: #negative
                   ReverseList[pileupcolumn.pos] += pileupread.alignment.query_sequence[pileupread.query_position]
             else:
                   ForwardList[pileupcolumn.pos] += pileupread.alignment.query_sequence[pileupread.query_position]

Note the 6th line, which implements a filter with a threshold of 10.


What are "split reads" and "intron clusters?" I'm working through the example data set for LeafCutter and the documentation mentions "split reads":

  This will cluster together the introns fond [sic] in the junc files listed
  in test_juncfiles.txt, requiring 50 split reads supporting each
  cluster and allowing introns of up to 500kb.

Searching "split reads" on Google yielded few results, and the results I did find were hard for me understand,

  Chimeric reads occur when one sequencing read aligns to two distinct
  portions of the genome with little or no overlap. Chimeric reads are
  indicative of structural variation. Chimeric reads are also called
  split reads.

One sequencing read aligns to two distinct portions of the genome, but then with little or no overlap? Isn't that a self-contradiction?

I can more or less deduce what "intron clusters" are but I figured I might as well try to solicit a more formal definition. I know what introns are and I can figure what it means to cluster them together. I guess I just want to know what split reads and intron clusters mean in the context of working out the example data for LeafCutter.

split reads - These are read that have two or more alignments to the reference from unique region of the read.  In this example a 150bp read sequenced from RNA could have base 1-75 aligning to the 3' end of exon2 and bases 76-150 aligning to the 5' end of exon3.  This would be a split read because it have two alignments (exon2 and exon3) and those alignment are from unique (non overlapping) parts of the read, 1-75 and 76-150.

Chimeric reads are an example of split reads.  Many people use "chimeric reads" to describe a type of noise that occur during PCR or sequencing where the bases from two different molecules are reported as one read and the alignments of that read become split when it aligns back to the location of the two original molecules.

"intron clusters" - In this case I think that the LeafCutter software is describing stacks of split reads to cluster to support the detection of an intron.


Generating burrows-wheeler-transform and a suffix array of a DNA sequence with less memory Also posted on biostars

After indexing

bwa index -a bwtsw reference.fa

I got files. like .bwt file and .sa file. The naive way of generating these file is:

Create all rotations of T
Sort the rotations lexicographically
Concatenate the last character of each rotation to form

BWT(T)

I used this naive way in my algorithm. The problem is it took too much memory.

My question is: is there an effective memory way to generate .bwt file and .sa file ?

Any help will be appreciated.

No, there is no way to decompress .bwt or .sa files, because they are not compressed. After indexing you'll get these files, which some programs (bwa) use. They should be in the same folder as your reference fasta file. What are you trying to do with these files? To get proper help here, please explain why you need these files, and show how you made the files.


How to read .bwt file .sa file? In BWA project in bwt.c file I found bwt_dump_bwt method and bwt_dump_sa method. I want to use them to read .bwt file and .sa file. I wrote the following program:

int main(void)
{
  char *tmp, *prefix;
    bwt_t *bwt;
    prefix = bwa_idx_infer_prefix("hg19bwaidx");

    tmp = calloc(strlen(prefix) + 5, 1);
    strcat(strcpy(tmp, prefix), ".bwt");
    bwt = bwt_restore_bwt(tmp);

  return 0;
}

I got the following error:

Undefined symbols for architecture x86_64:
  "_bwa_idx_infer_prefix", referenced from:
      _main in readbwt-3022ec.o
  "_bwt_restore_bwt", referenced from:
      _main in readbwt-3022ec.o
ld: symbol(s) not found for architecture x86_64
clang: error: linker command failed with exit code 1 (use -v to see invocation)

Any help ?

Thanks.

You need to include bwa.o in your call to clang, since it holds the definitions of those functions. This will cause a good number of object files, so you might as well just clang -o foo *.o foo.c -lz -lpthread -lm -lrt and get it over with. Since this was quite obviously never meant to be done, you'll have to remove at least one main() function from bwa.c.

Note that further questions regarding this will be considered off-topic and closed, since it's purely about compiling something that you want to write yourself.


ANNOVAR, 1000g database calling I have recently downloaded Annovar, and succesfully filtered my files with a few databases, downloaded from Annovar. But when I try to do the exact same thing with the 1000 genomes project database, I get this error many times in the terminal, one for each line of **.

Argument "&" isn't numeric in numeric eq (==) at annotate_variation.pl line 1962, <DB> line &&.

& can be A, C, G, or T.

&& every line of the database

I have tried filtering with 1000g, 1000g2010 and 1000g2010jul databases, but all give these errors.

An example of the commands I use for filtering:

perl annotate_variation.pl --filter --dbtype CEU.sites.2010_07 testfile1.txt humandb/

Can anyone tell me if I'm doing something wrong, if this is known, and if I can do anything about it?

That's a warning. Presumably, a match failed ($& is the default Perl variable which holds the result of a match operation). So something in your input is not exactly what the script was expecting. You can probably safely ignore these messages.

Alternatively, update your question and add a few example lines of input which we can use to reproduce the error and maybe we can figure something out. If that isn't possible, at least show us line 1962 of the annotate_variation.pl script since that's the one that's failing.


How can I use my Myheritage DNA results file for further analysis? I had my DNA tested by Myheritage and they sent me a csv file with RSID, Chromosome, position and result (which base) with about 700,000 rows.
I understand most analyses of DNA use VCF files but is there anything i can do with this csv file i.e. check for genetic health-related genes?
I am not a bioinformaticist but I am a scientist and I can use R and python.  Ive heard of the Gnomad database but not sure if i can match things to my csv file.

A csv file with RSID, chromosome, position and result is enough for what you want to do and these are the core columns of a VCF (which is just a TSV with some headers describing how it was made).
Given you have ~700k rows I suspect that yes there will be genetic health-related SNVs (single nucleotide variant). Disclaimer; these data could include info pertaining to your health and your family's health, I strongly advise that you speak to a genetic counsellor to understand these things.
Gnomad is a population database so generally speaking (plenty of exceptions of course) if there is a SNV in there more than a few times its probably not causing disease. https://www.ncbi.nlm.nih.gov/clinvar/ and https://omim.org/ are example of disease databases.
Here's a pandas example with clinvar (https://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/ you'll need to match your version).
import pandas as pd
df = pd.read_csv(
    'clinvar.vcf.gz',
    sep='\t',
    comment = '#',
    header=None,
    names = ['CHROM','POS','ID','REF','ALT','QUAL','FILTER','INFO']
)
df.head()

CHROM   POS ID  REF ALT QUAL    FILTER  INFO
0   1   930188  846933  G   A   .   .   ALLELEID=824438;CLNDISDB=MedGen:CN517202;CLNDN...
1   1   930203  972363  C   T   .   .   ALLELEID=959431;CLNDISDB=MedGen:CN517202;CLNDN...
2   1   930248  789256  G   A   .   .   AF_ESP=0.00347;AF_EXAC=0.00622;AF_TGP=0.00280;...
3   1   930275  969662  T   G   .   .   ALLELEID=959432;CLNDISDB=MedGen:CN517202;CLNDN...
4   1   930336  843786  G   A   .   .   ALLELEID=824439;CLNDISDB=MedGen:CN517202;CLNDN...

Then lets say one of your 700k SNVs is at chrom 1, position 930336 and you have an A
print(*df.query('(CHROM == 1) & (POS == 930336) & (ALT == "A")')['INFO'].str.split(';'))

Gives a list from which you could pull out CLNSIG which here is Uncertain_significance
['ALLELEID=824439', 'CLNDISDB=MedGen:CN517202', 'CLNDN=not_provided', 'CLNHGVS=NC_000001.11:g.930336G>A', 'CLNREVSTAT=criteria_provided,_single_submitter', 'CLNSIG=Uncertain_significance', 'CLNVC=single_nucleotide_variant', 'CLNVCSO=SO:0001483', 'GENEINFO=SAMD11:148398', 'MC=SO:0001583|missense_variant', 'ORIGIN=1']

So you could easily parse your 700k SNVs and get their respective clinical info this way. Please think about this carefully before you do it though, there are many implications including life insurance and mental health! I'm of the opinion that this sort of data is useful but MANY people disagree. Further, this advice is general in nature only, I accept no liabilty for what you or anyone else chooses to do with this publicly available data.


Can anyone point me towards a HiC Dataset with domains of high confidence? is there a dataset online containing a HiC contact matrix and file with TAD boundary locations, which were generated, either by hand or automated, with high accuracy? Although I already used programs such as SpectralTAD or GMAP to generate TADs myself, a dataset which was computed with higher accuracy would be nice to check my project against.
Preferred attributes would be:

human genome
whole genome
10 kb resolution or higher

Thank you very much in advance!
edit: One dataset I found is http://chromosome.sdsc.edu/mouse/hi-c/download.html, which includes multiple 10 kb datasets with TADs called.

I guess your best chance is one of the in-situ HiC papers from the Lieberman-Aiden lab. There are not too many HiC datasets with 10kb-or-better resolution out there for human. I would check this one:
A 3D Map of the Human Genome at Kilobase Resolution Reveals Principles of Chromatin Looping


How to concatename file names in Snakemake? I want to pass a series of files (e.g., a, b) to a program prog like: prog a,b.

Is there a way to do this in Snakemake?

My attempt is:

rule run_prog:
  input: "a", "b"
  output: "output"
  run: "prog {input}"

But the actual command Snakemake generated is prog a b instead of prog a,b.

Any solution to pass a comma concatenated file names to 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)}"


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

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

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

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

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

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

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

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

I try the below code, which did not work.

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

Thanks in advance.

Ok I have figured it out and I am sharing this in case anyone can get help. I used
sce <-addPerCellQC(sce)

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


How to tell if library generation for RNA-seq experiments are stranded vs not stranded I would like to use RNA-seq data from the NIH GEO database. How can I tell if library generation for the RNA-seq experiments are stranded vs not stranded?

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


What are the count units in DEseq2? This is a pretty silly/simple question, I suppose.

What units are DESeq2 counts? Or are the units arbitrary but internally consistent estimate of actual reads?

The counts are "reads" for single-end datasets and "fragments" for paired-end datasets. In other words, they're what featureCounts and htseq-count produce. The "normalized counts" that you'll be able to access in DESeq2 are, as aptly named, simply normalized versions such that values are comparable across samples.


Detecting SARS-CoV-2 variants from the mixed virus population I have a fastq file from  Nanopore sequencing data that contains reads from the UK and South Africa variants of SARS-CoV-2. These variants are identified by three key mutations in the receptor-binding domain. I want to determine the number of reads then belongs to UK and South Africa variant. I will create a bam file from the fastq file using the SARS-CoV-2 reference genome. Then, how do I determine the reads belonging to each variants? Are there any tools available to classify the reads based on the pre-defined mutation criteria?

Based on a quick request from someone else, I wrote some code a couple of years ago for assigning read groups based on variants at particular locations; this was set up for HIV variant splitting, so I expect it should also work on SARS-CoV-2 nanopore reads as well:
https://gitlab.com/gringer/bioinfscripts/-/blob/master/samVarSplitter.pl
It processes a SAM stream [i.e. in text mode] using Perl, so might be a bit slow, but should be fine for a few hundred thousand sequences. You'd use it as follows:
samtools view -h mapped_reads.bam | \
  samVarSplitter.pl [-ref <refName>] [-pos <int>] | \
  samtools sort > grouped_reads.bam

Then the reads can be split into different files based on their group using samtools split:
Usage: samtools split [-u <unaccounted.bam>] [-h <unaccounted_header.sam>]
                      [-f <format_string>] [-v] <merged.bam>


CDS length for each human gene Does anyone know where and how could I download a list of all human genes and the length of the coding sequence for each gene? Is it possible to do this on the NCBI site, ensembl?

While I haven't found a way to limit the results to the canonical transcript only, you can get a list of genes, transcripts and their CDS lengths using Ensemble's BioMart. I have already set it up for you, you can see the results, and modify them, here (click on the "Results" link if you don't see them).

Essentially, you just need to go to BioMart, and

select "Ensembl Genes 96" (the number will change if the version changes) as the database and "uman Genes" as the dataset.
Click on "Filters", and set Gene type to coding and Transcript type to protein_coding.
From "Attributes", select whatever you want to see. The "CDS Length" is under "Structures".


Modifying a SAM header after removing all non-primary reads I subset a BAM to only include primary reads using the following samtools commands:

samtools view -F 256 input.bam > input.primaryOnly.sam

Now, in order to convert this SAM into a new subsetted BAM, I would need to add a header. The way I normally do this is with Picard:

java -jar picard.jar ReplaceSamHeader \
      I=input.primaryOnly.sam \
      HEADER=input.header.sam \
      O=new.sam

whereby the header input.header.sam is from the original BAM:

samtools view -H input.bam > input.header.sam

Then I could convert the SAM into a BAM with

samtools view -S -b i new.sam > new.bam

My problem is, I'm using the original BAM header. If I look back at the BAM, it now contains wrong information, e.g.

samtools flagstat new.bam

might save that secondary alignments exist, when they don't.

How do I get a BAM header that makes sense given my new BAM (with only primary alignments)?

Your first command only needs a slight modification to add in -h. This will create a SAM file with a header.

samtools view -h -F 256 input.bam > headered.input.primaryOnly.sam

  -h        Include the header in the output.

If your primary end goal is to create a BAM, then as Pierre has pointed out you can create BAM files directly from samtools view with the -b option:

samtools view -b -F 256 input.bam > input.primaryOnly.bam

And if your end-goal is a sorted BAM file, then pipe the BAM output through samtools sort:

samtools view -b -F 256 input.bam | samtools sort > sorted.input.primaryOnly.bam

So that's fixing what I think your problem is, but what about the specific question that you have asked, i.e. how can a header be added to an existing headerless SAM file?

If you have an existing BAM file, you can use Picard (as you have mentioned), or manually piece the files together with cat, which will have a similar function (although Picard probably includes a few checks to make sure that the header is consistent with the input file):

(samtools view -H input.bam; cat input.primaryOnly.sam) > headered.input.primaryOnly.sam

If the only thing that matters is that the BAM file has correct reference information (which is fairly common), then the -T option can be used with samtools view. Again, -h or -b is required to make sure both the header and the usual SAM fields are displayed:

samtools view -T reference.fa -h input.primaryOnly.sam > headered.input.primaryOnly.sam
samtools view -T reference.fa -b input.primaryOnly.sam > headered.input.primaryOnly.bam


How could I match atom orders between a .mol2 and a .pdb? When I was drawing conformation of Autodock result, I noticed that Autodock generated pdbqt from dlg has some structural problem (e.g. Benzene ring missing), and wants to correct this. However I noticed that the atom order of pdbqt doesn't match with the original mol2's, and caused some problems when I tried to fix.

So how could I reorder pdbqt according to mol2's atom coordinate?

Atom reordering is a common problem in compchemistry.

Rdkit (a python package) can do this, but it is limited by the formats it can read and mol2 files are a bit hit or miss. It works really well with SMILES, SMARTS and mol (sdf) files. But the writing may cause problems with Brookhaven pdb and mol2 files.
So the formats you have are both problematic. One is a pdbqt (a PDB file with partial charges (Q, generally Gasteiger calculated) and atom type (T) and the other a mol2. Let's pretend first that they are not.

Also the hydrogens might be made implicit depending on your mol2 file.

from rdkit import Chem
ref = Chem.MolFromMol2File('ligand.mol2')
target = Chem.MolFromPDBFile('ligand.pdb')
atoms = target.GetSubstructMatches(ref) #: List[int]
new = Chem.RenumberAtoms(target, atoms)
Chem.MolToMolFile(new, 'neoligand.mol')

Extra steps would be required for the protons and if your starting ligand was not from a file, but a SMILES.

mol = Chem.AddHs(mol) #protonate explicitly
Chem.GetSSSR(mol) #resonance fixing
AllChem.EmbedMolecule(mol) #initialise for 3d.
AllChem.UFFOptimizeMolecule(mol, maxIters=2000)

So this is all nice, but the formats may not work. Luckily openbabel can deal with this problem.

obabel -i mol2 input.mol2 -o pdbqt -O output.pdbqt

However, one thing to keep an eye out is what will be recalculated. Say, sdf files do not have partial charges or what PDB calls "atom names" and Rdkit calls "atom labels" (About the latter, I actually wrote a script to relabel atom labels in pdbs from mol files), while mol2 and PDBQT have.


Genetic Genealogy related to SNP variants and INDEL variants Background:
I am looking for some advice on the use of SNP variants versus INDEL variants for determining Genetic Genealogy, specifically for determining accurately the "START" and "END" of shared DNA segments (shared between two DNA testers' results). The test data used is RAW DNA from direct to consumer DNA testing companies such as AncestryDNA and 23andMe.

Question:
For the purpose of Genetic Genealogy, are the INDEL variants as useful as the SNP variants? Should the INDEL variants be ignored? Or, alternatively, do the INDEL variants hold a greater significance than the SNP variants?

I hope that using INDEL variants would increase the data in the comparison algorithm, thus increasing the certainty and resolution when determining each shared segment.

InDels are just as informative, but are less reliable to call. Errors should be averaged out in this kind of analysis though (so long as they don't predominate).

Both SNPs and InDels are valid data points to consider when comparing triads, but identifying IBD segments will only ever be probabilistic.


How to represent trans-spliced genes in GTF? For example, see this gene (nad1) in ENA:
http://www.ebi.ac.uk/ena/data/view/ABI60879

If you look at the XML for that gene you see the following:

join(
             DQ984518.1: 324706 .. 325091 ,
  complement(DQ984518.1:  24417 ..  24498),
  complement(DQ984518.1:  22828 ..  23019),
             DQ984518.1:   3484 ..   3542 ,
  complement(DQ984518.1: 153702 .. 153960)
)

Which shows 5 exons joined out of phase and out of order. Is there a valid GTF representation of this?

How to dump a 'non-canonically spliced' gene into GTF? i.e. what's the recommendation?

Cross-posted on biostars

To my knowledge there's no defined way to deal with that in GTF. GFF3 handles trans-splicing (you'll have to scroll down to "trans-spliced transcript") by giving an individual transcript multiple parents (e.g., ID=some_transspliced_gene;Parent=gene1,gene2). You could use the same methodology with GTF files, but just note that it'll break most downstream programs.


Human body sites knowledge base I'm looking for a lightweight knowledgebase describing the human body to annotate disease sites.
I do not need a great level of detail, I just need a kind of basic organs/sub-organs taxonomy. I checked out the Foundational Model of Anatomy ontology, but it's over-kill for my purposes and really cumbersome (Protege fails when it tries to load FMA.owl).

Is there a simple online knowlegebase available for this?

I am the developer of Uberon and I would be happy to help you with what you need, either from Uberon, or from the FMA.

You mentioned you need something simpler than FMA. There are a variety of tools for creating custom subsets. Additionally, some ontologies provide ready-made subsets for particular purposes.

For example, we make a subset called 'basic.obo' from Uberon:
http://uberon.github.io/downloads.html#subsets

You can examine this directly in either oboedit or Protege5 (note Protege5 parses obo so you don't need to convert).

From your question, you didn't indicate you needed programmatic access, but if you do there are a variety of options, as well as more programmer-friendly JSON exports.

Most ontologies have a number of different options for web-browsing, you can see some of the options for Uberon here:
http://obofoundry.org/ontology/uberon.html

You mentioned annotating disease sites. There may be a tools (either desktop or web-based) for helping with this - e.g. providing autocomplete over a desired subset of an ontology.

You may also want to ask your question on a mailing list like obo-discuss. Unfortunately many ontology developers don't read stackexchange, but most are always keen to hear of community requirements, to help them provide better and more useful products.


How drastic would the result be if I use fasta or reference assembly from ucsc and gtf from gencode? There are difference annotation file for UCSC and gencode.

But if I use the reference assembly from UCSC and the GTF from Genocode or vice versa would my downstream results would be wrong?

Never use genomes or annotations from UCSC, they're poorly versioned and only recently with mouse and human have they even included all of the contigs. For fasta/GTF files from early in the GRCh38 release, you can tell whether you're using UCSC or Gencode by the presence/absence of _random contigs, which will only exist for UCSC. These were mostly later split into the actual contigs, so recent download from UCSC should more closely match what you find at Gencode/Ensembl. Further, that time predated UCSC beginning to adopt Gencode's vastly superior annotations, so if your GTF file has instances where the same gene ID is on either multiple strands or multiple chromosomes (this is obviously biologically impossible) then you have a UCSC GTF file.

In general, with early GRCh38 releases your only real issues will be with _random contigs, which are a minority of the genome and don't have all that many genes. But really, you should be keeping track of the sources of your files and ensuring that they're compatible.

Update: I should expand a bit on my "Never use genomes or annotations from UCSC" comment. In point of fact the genomes themselves aren't so terrible. Early on UCSC had the bad habit of concatenating contigs together into _random "chromosomes", but they seem to have mostly kicked that habit as of late. Note, however, that there are no versions for their genomes. Since reference genomes continue to get updates over time (mostly through the addition of patches) the lack of actual versions means you have to manually check if a recently downloaded file matches what may have been downloaded either previously or by someone else. This has obvious consequences when it comes to reproducibility. The same issue occurs for annotations from UCSC, but they have the additional problem of historically having biologically incoherent concepts of genes. That is, they will contain the same gene in multiple places with multiple orientations, which will break many tools in both obvious and completely unclear ways. For example, DEXSeq will simply break with an error message if given a UCSC annotation, since they break biological plausibility. If you were to use these annotations files with deepTools, you wouldn't get an error message, but the resulting output would be only partial, due to the biologically impossible annotation effectively corrupting most obvious ways of storing annotation data in a data structure (i.e., you can no longer treat IDs as unique). This could have downstream ramifications on biological interpretation of results.


Calculate the percentage of each unique phylogenetic tree in a BEAST output I have a nexus formatted BEAST output containing 20,000 phylogenetic trees of seven taxa. Is there any way to get the percentage of each unique phylogenetic tree contained in this output?

I already made an unsuccessful attempt with R.

I finally managed to do it in R. Here is my code:

install.packages('devtools')
library(devtools)
install_github('santiagosnchez/rBt')
library(rBt)

beast_output <- read.annot.beast('beast_output.trees')
beast_output_rooted <- root.multiPhylo(beast_output, c('taxon_A', 'taxon_B'))
unique_topologies <- unique.multiPhylo(beast_output_rooted)

count <- function(item, list) {
  total = 0
  for (i in 1:length(list)) {
    if (all.equal.phylo(item, list[[i]], use.edge.length = FALSE)) {
      total = total + 1
    }
  }
  return(total)
}

result <- data.frame(unique_topology = rep(0, length(unique_topologies)),
                     count = rep(0, length(unique_topologies)))
for (i in 1:length(unique_topologies)) {
  result[i, ] <- c(i, count(unique_topologies[[i]], beast_output_rooted))
}

result$percentage <- ((result$count/length(beast_output_rooted))*100)


Cellranger gives error I am trying to run cellranger but I get fastq permission denied error
fi1d18@RBGO-Server2:/opt$ cellranger-atac count --id test --reference refdata-cellranger-arc-GRCh38-2020-A-2.0.0 --fastqs /data/fi1d18/Downloads/FASTQ_scATAC/ --chemistry ARC-v1
error: couldn't open MRO file
    Caused by: Permission denied (os error 13)
fi1d18@RBGO-Server2:/opt$

I have added the path by
export PATH=/opt/cellranger-atac-2.0.0:$PATH

And by vim I have also added the path
fi1d18@RBGO-Server2:/opt$ vim ~/.bashrc

Any help please?

Two of those command line parameters should be accessible locations (see here). Could you please add the output of the following commands:
ls -l /data/fi1d18/Downloads/FASTQ_scATAC/ | head
and
ls -l refdata-cellranger-arc-GRCh38-2020-A-2.0.0 | head
I suspect that the refdata link should be an absolute path, rather than a relative path, i.e. something like this:
cellranger-atac count --id test --reference /opt/refdata-cellranger-arc-GRCh38-2020-A-2.0.0 --fastqs /data/fi1d18/Downloads/FASTQ_scATAC/ --chemistry ARC-v1


Attractor Landscape Analysis I have come across a modeling toolbox, ATLANTIS, which is able to determine cell fates in silico based on the input models provided. This MATLAB-based toolbox is built on a method called "Attractor Landscape Analysis." I tried to find a paper dedicated to explaining this method but I did not have any luck so far. Can any of you folks please explain it to me, especially the term "attractor"?

From Wikipedia:
In the mathematical field of dynamical systems, an attractor is a set of numerical values toward which a system tends to evolve, for a wide variety of starting conditions of the system.

In biology you can think of an attractor as a state towards the system is moving to. A good example is the process of differentiation where stem cells compromise to more specialized cells, that could be thought as the attractors of the process.


How to version the code and the data during the analysis? I am currently looking for a system which will allow me to version both the code and the data in my research.

I think my way of analyzing data is not uncommon, and this will be useful for many people doing bioinformatics and aiming for the reproducibility.

Here are the requrements:

Analysis is performed on multiple machines (local, cluster, server).
All the code is transparently synchronized between the machines.
Source code versioning.
Generated data versioning.
Support for large number of small generated files (>10k). These also could be deleted.
Support for large files (>1Gb). At some point old generated files can  permanently deleted. It would be insane to have transparent synchronization of those, but being able to synchronize them on demand would be nice.

So far I am using git + rsync/scp. But there are several downsides to it.

Synchronization between multiple machines is a bit tedious, i.e. you have to git pull before you start working and git push after each update. I can live with that.
You are not supposed to store large generated data files or large number of files inside your repository.
Therefore I have to synchronize data files manually using rsync, which is error prone.

There is something called git annex. It seems really close to what I need. But:

A bit more work than git, but that's ok.
Unfortunately it seems it does not work well with the large number of files. Often I have more that 10k small files in my analysis. There are some tricks to improve indexing, but it doesn't solve the issue. What I need is one symlink representing the full contents of directory.

One potential solution is to use Dropbox or something similar (like syncthing) in combination with git. But the downside is there will be no connection between the source code version and the data version.

Is there any versioning system for the code and the data meeting the requirements you can recommend?

There is a couple of points to consider here, which I outline below. The goal here should be to find a workflow that is minimally intrusive on top of already using git.

As of yet, there is no ideal workflow that covers all use cases, but what I outline below is the closest I could come to it.

Reproducibility is not just keeping all your data

You have got your raw data that you start your project with.

All other data in your project directory should never just "be there", but have some record of where it comes from. Data processing scripts are great for this, because they already document how you went from your raw to your analytical data, and then the files needed for your analyses.

And those scripts can be versioned, with an appropriate single entry point of processing (e.g. a Makefile that describes how to run your scripts).

This way, the state of all your project files is defined by the raw data, and the version of your processing scripts (and versions of external software, but that's a whole different kind of problem).

What data/code should and should not be versioned

Just as you would not version generated code files, you should not want to version 10k intermediary data files that you produced when performing your analyses. The data that should be versioned is your raw data (at the start of your pipeline), not automatically generated files.

You might want to take snapshots of your project directory, but not keep every version of every file ever produced. This already cuts down your problem by a fair margin.

Approach 1: Actual versioning of data

For your raw or analytical data, Git LFS (and alternatively Git Annex, that you already mention) is designed to solve exactly this problem: add tracking information of files in your Git tree, but do not store the content of those files in the repository (because otherwise it would add the size of a non-diffable file with every change you make).

For your intermediate files, you do the same as you would do with intermediate code files: add them to your .gitignore and do not version them.

This begs a couple of considerations:

Git LFS is a paid service from Github (the free tier is limited to 1 GB of storage/bandwidth per month, which is very little), and it is more expensive than other comparable cloud storage solutions. You could consider paying for the storage at Github or running your own LFS server (there is a reference implementation, but I assume this would still be a substantial effort)
Git Annex is free, but it replaces files by links and hence changes time stamps, which is a problem for e.g. GNU Make based workflows (major drawback for me). Also, fetching of files needs to be done manually or via a commit hook

Approach 2: Versioning code only, syncing data

If your analytical data stays the same for most of your analyses, so the actual need to version it (as opposed to back up and document data provenance, which is essential) may be limited.

The key to get this this working is to put all data files in your .gitignore and ignore all your code files in rsync, with a script in your project root (extensions and directories are an example only):

#!/bin/bash
cd $(dirname $0)
rsync -auvr \
    --exclude "*.r" \
    --include "*.RData" \
    --exclude "dir with huge files that you don't need locally" \
    yourhost:/your/project/path/* .

The advantage here is that you don't need to remember the rsync command you are running. The script itself goes into version control.

This is especially useful if you do your heavy processing on a computing cluster but want to make plots from your result files on your local machine. I argue that you generally don't need bidirectional sync.


RNA velocity: competing explanations for variable ratio of spliced to unspliced transcript In "RNA velocity of single cells", La Manno et al. look at ratios between spliced and unspliced mRNA as a way of estimating the velocity of a cell through transcriptome space.

In checking for a universal steady-state spliced-unspliced ratio, they write "We found that the steady-state behaviour of most genes across a wide range of cell types was consistent with a single fixed [ratio between spliced and unspliced] (Extended Data Fig. 3a–c). However, 11% of genes showed distinct slopes in different subsets of tissues (Extended Data Fig. 3d, e), suggesting tissue-specific alternative splicing (Extended Data Fig. 3f) or degradation rates."

Tissue-specific alternative splicing or degradation rates certainly could account for these patterns, but how do La Manno et al. rule out that some tissues are simply not in a steady state? The figure they cite (screenshotted below) uses data from the Tabula Muris adult mouse cell encyclopedia. In an adult mouse, not all cells would not be expected to be in steady state -- for example, T-cell progenitors from the thymus and epithelial cells in many organs undergo continual replenishment.

If cells aren't in a steady state then one wouldn't expect to see the results in (e), where there are two slopes. One would instead expect to see results more like Plekha3 in panel (c), where there's more variability around the the steady state line (or perhaps no steady state line at all).


How low does an E-Value need to be for a gene to code for the same protein as query? I'm looking at a bacterial genome and trying to find specific PEGs (Protein Encoding Genes) by BLAST searching. I'm taking the genes found with the lowest E-value in order to psiBLAST them to check it codes for the same protein. I'm not exactly sure of how high the E-Value needs to be for me to conclude that the protein is not coded for by the subject gene.
Does anyone have any advice on how to make this process a little more efficient/scientific?
Thank you.
(By the way, I'm using BLAST with a BLOUSUM62 matrix.)

I assume that you mean "is homologous to" by "encodes the same protein as". Identical ("same") proteins should be fairly trivial to identify using a tool such as tblastx.
You should probably not use E-value for this. E-value is similar to a p-value, in that it represents some aspect of the data (in this case a sequence and a database) compared to some null hypothesis (meaning in this case statistics on the database). Therefore, it depends both on the specifics of the data and the hypothesis, and there is no universal threshold that will work across all sequences and all databases. (See discussion of the meaning of e-value.)
E-value does not in any way represent the content of the sequence or its relationship to a sequence in the database. You could have E-values close to zero or E-values > 10 for a perfect sequence match or a very loose sequence match, depending on the input sequence and the database.
However, BLAST output does include other information that might be helpful. You could (for example) look at the length of the alignment and the number of mismatches/identities in the alignment. This will give you a pretty good guess at how similar the two sequences are; unless they are very divergent you would expect them to align across most of the length with somewhat high protein identity or similarity (>25% probably?). This process will still not work 100% of the time, but I think it is closer to your intent than any given approach involving E-values.


using snakemake shadow rules to store temp files on local nodes Snakemake is not really designed to store temporary files to local disk, but in the documentation, they suggest a workaround using "shadow rules".
The problem I have is, that our front node does not have /scratch/.. directory and therefore snakemake complains right away
snakemake my_result.tsv -p --jobs 10 --shadow-prefix /scratch/$USER --cluster "qsub "
   ...
   PermissionError: [Errno 13] Permission denied: '/scratch'

Every single computational node does have the directory, only the frontend does not. Is there any easy workaround?
--- edit ---
I will try to clarify where I think is the problem.
The head node does not have anything mounted to /scratch, I think that's why I am getting that error. I sort of expected that snakemake will attempt do create the shadow environments for individual executed jobs, because if it does, everything would be fine. All the computing nodes do have their local disks mounted to /scratch.
Of course, I could use a different shadow directory, but all the directories that are mounted both at the front-end and computing nodes are NOT local disks of individual nodes, which would cause more network traveling of the temporary data.
--- edit 2 ---
Thanks @Chris_Rands for finding the relevant lines. I opened an Issue on Snakemake GitHub. Once I will know more, I will get back with the answer here.

What is happening is your Admin is not allowing you access to /scratch ... which is undertandable. Essentially you are attempting to access the headnode /scratch. If every user on the cluster had access, it would be a bit chaotic.
I am surprised by the error, not that I would use Snakemake I'd code it, because you would expect a data pipeline tool as established as Snakemake to have a solution to working on a cluster. What you are describing is a standard cluster set-up between headnode and 'sub-node' (the technical term is slavenode, but I hate that term). The headnode is mirrored across the cluster, hence the error. Once you submit via qsub the mirror will 'display' the headnode layout on the subnode.

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

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


Error in annotating row names as data points in a PCA plot I made a PCA plot and was trying to plot the eigenvenctors in R so that each data point is actually the sample name. All the samples are in column 1. The following are the R commands I am using:

    fn <- "Plot.txt"
    evecdat <- read.table(fn, col.names=c("Sample", "PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8", "PC9", "PC10", "pc11","pc12", "pc13", "pc14", "pc15","pc16", "pc17", "pc18", "pc19", "pc20", "pop", "color"))
    jpeg("Plot1.jpeg", width = 20, height = 20, units = 'in', res = 300)
    identify(evecdat$PC1, evecdat$PC2  = NULL, labels = seq_along(evecdat$Sample), pos = TRUE, n = length(evecdat$Sample), plot = FALSE, atpen = FALSE, offset = 0.5, tolerance = 0.25, …)

    plot(evecdat$PC1, evecdat$PC2, col=as.vector(evecdat$color), xlab="PC1", ylab="PC2", pch="24", cex.lab=2.5, cex=1.5, cex.axis=1)
    legend("topright", ncol=1, legend=c("X", "Y", "Z"), col=c("Blue", "Red", "Green"), pch="24", cex=2, bty = "o")
    while (!is.null(dev.list())) dev.off()

Here, all the row names are under the Sample column name that I defined earlier.

Maybe try

plot(evecdat$PC1, evecdat$PC2, xlab="PC1", ylab="PC2", pch="24", cex.lab=2.5, cex.axis=1, type="n")
text(evecdat$PC1, evecdat$PC2, evecdat$Sample, col=as.vector(evecdat$color))

This should plot the sample names as text (taken from the evecdat$Sample column), colored as specified in the evecdat$color column. See ?text for details.


Path to career in bioinformatics? I've been completing various courses on bioinformatics through edX and Coursera, and I've made a strong footing in the material.

I'm not sure what to pursue next in bioinformatics that would bring me closer to research in the field or employment.

I don't have bioinformatics at my university and thus am unsure of how to progress.

I have a background in mathematics, chemistry, physics and computer science. So, what would you recommend I do next?

The Bioinformatics Stronghold at Rosalind is an excellent resource for building your understanding of many concepts and algorithms you'll see frequently in bioinformatics and genomics.

I'd also say that, if you can make it work, an internship where you can "get your hands dirty" with real problems is a very valuable experience. This option may not be available to everyone due to a variety of factors, but when it is an option it provides a tremendous opportunity to build experience with little risk to the intern or the employer.


Illustrate a 3D visualisation of the three main PCs using plot3d() package in R? I have a dataset (view/download) here. I need to run PCA on this data set and need to illustrate a 3D visualisation of the three main PCs using plot3d() package.

I have looked for more information before coming here on which package to use, some people use princomp and some us prcomp (in university we have personally been using prcomp). The confusing part is were to move forward with this?

The full question is "Illustrate a 3D visualisation of the three main PCs using plot3d() package in R. Use the following colours for samples (data points) belonging to the various immune subgroups/subtypes (C1: red, C2: yellow, C3: green, C4: cyan, C5: blue, C6: purple)"

After reading in the CSV and loading rgl library

In the current state I have removed the rows with the subgroups: see code below

# Remove last row from dataset
dataMinusLastRow = data[1:(nrow(data)-1),]

# Convert dataset to numeric and apply
coltemp <- colnames(dataMinusLastRow)
pcaData = t(apply(dataMinusLastRow, 1, as.numeric))
colnames(pcaData) <- coltemp

# PCA
pc = prcomp(t(pcaData), scale=T, center=T)
summary(pc) # From this we can see that PC184 accounts for 95% of cumulative proportion

# Getting the cumulative proportion in numeric dataset
vars <- apply(pc$x, 2, var)
props <- vars / sum(vars)
pc.cs <- cumsum(props) # Cumulative Proportion
pc.index<-min(which(pc.cs>0.95)) # Index of which cumulative proportion > 95%

I think the main part which confuses me is that it asks for the three main PC's which I am not sure what they are and have had not much support from university with information regarding this.

Thanks to anyone who can help in any shape/form!

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


How to simulate NGS reads, controlling sequence coverage? I have a FASTA file with 100+ sequences like this:

>Sequence1
GTGCCTATTGCTACTAAAA ...
>Sequence2
GCAATGCAAGGAAGTGATGGCGGAAATAGCGTTA
......

I also have a text file like this:

Sequence1 40
Sequence2 30
......

I would like to simulate next-generation paired-end reads for all the sequences in my FASTA file. For Sequence1, I would like to simulate at 40x coverage. For Sequence2, I would like to simulate at 30x coverage. In other words, I want to control my sequence coverage for each sequence in my simulation.

Q: What is the simplest way to do that? Any software I should use? Bioconductor?

I am working on a Illumina sequencing simulator for metagenomics: InSilicoSeq

It is still in alpha release and very experimental, but given a multi-fasta and an abundance file, it will generate reads from your input genomes with different coverages.

From the documentation:

iss generate --genomes genomes.fasta --abundance abundance_file.txt \
    --model_file HiSeq2500 --output HiSeq_reads

Where:

# multi-fasta file
>genome_A
ATGC...
>genome_B
CCGT...
...

# abundance file (total abundance must be 1!)
genome_A    0.2
genome_B    0.4
...

I didn't design it to work with coverage but rather abundance of the genome in a metagenome, so you might have to do a tiny bit of math ;)


How to query the Human Microbiome Project (HMP) to find all subjects with both 16s and WGS workups? I am looking for a query to run on the HMP database that will return all subjects who have had BOTH 16s and whole genomes sequence (WGS) workups.  I am currently using this query...

file.matrix_type = 16s_community OR file.matrix_type = wgs_community AND (file.format = FASTA OR file.format = FASTQ)

and sifting through the results by adding them to my shopping cart, browsing the file types until I have added all of the samples from one subject, and then I check to see if I have at least one 16s and one WGS.  So far, I have not found one person who has had both 16s and WGS where both files were publicly accessible, and most do not have WGS at all.  There must be a way to get subjects with both workups done to make this faster.

After trying a few different queries, I don't believe there is a more specific one than this and I found that by searching and clicking through all of the results given by this query that eventually I found patients that fit the criteria of having both WGS and 16s.


Where can I find .cdf files for Affymetrix ChIP-chip data analysis? I am analysing microarray ChIP data generated several years ago. I want to use the affy package for R. However, in order to load the CEL files and proceed with the analysis, I need the pombe annotation file which comes under the name of Sp20b_M_v04.cdf. I have looked everywhere for this file but I am unable to locate it. Is anyone in possession of this file or know where I can download it from?

It looks like it's in GEO here: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL10187

There is apparently a zipped copy linked at the bottom in the supplementary files.

When people upload a study based on an array, there is a corresponding 'platform' record where the array details should be deposited.

Edit: you will probably need to gunzip the file first, to get it into CDF format.


How to convert VCF format file into Fasta(.fa) with the Reference data I need a Software or Python Program for Converting VCF(.vcf) to  FASTA(.fa) format with the help of reference (.fa)

I used FastaAlternateReferenceMaker previously. You will need to download GATK first. Next you prepare your reference genome (reference.fasta) and your vcf file (input.vcf). You will call the GATK function like this:

 java -jar GenomeAnalysisTK.jar \
   -T FastaAlternateReferenceMaker \
   -R reference.fasta \
   -o output.fasta \
   -V input.vcf \

output.fasta will contain the new fasta with snps inserted at sites specific by the vcf file


How to remove duplicates from a fasta file using python I am using the following command:
from Bio import SeqIO
import sys
import re

    fasta_file = (sys.argv[1])
    for myfile in SeqIO.parse(fasta_file, "fasta"):
      if len(myfile) > 250:
       gene_id = myfile.id
       mylist = re.match(r"H149xcV_[^\W_]+_[^\W_]+_[^\W_])_[^\W_]+", gene_id)
       print (">"+list.group(1))

And this is providing me with duplicates of the same gene:
>H149xcV_Fge342_r3_h2
>H149xcV_bTr423_r3_h2
>H149xcV_kN893_r3_h2
>H149xcV_DNp021_r3_h2
>H149xcV_JEP3324_r3_h2
>H149xcV_JEP3324_r3_h2
>H149xcV_JEP3324_r3_h2
>H149xcV_SRt424234_r3_h2
>H149xcV_Fge342_r3_h2
>H149xcV_Fge342_r3_h2

How can I reformat my command so that I only receive unique gene id's:
>H149xcV_Fge342_r3_h2
>H149xcV_bTr423_r3_h2
>H149xcV_kN893_r3_h2
>H149xcV_DNp021_r3_h2
>H149xcV_JEP3324_r3_h2
>H149xcV_SRt424234_r3_h2

from Bio import SeqIO
import sys
import re

unique = []
fasta_file = (sys.argv[1])
for myfile in SeqIO.parse(fasta_file, "fasta"):
    if len(myfile) > 250:
        gene_id = myfile.id
        if gene.id not in unique:
            unique.append(gene.id)
            mylist = re.match(r"H149xcV_[^\W_]+_[^\W_]+_[^\W_])_[^\W_]+", gene_id)
            print (">"+mylist.group(1))

```


In SBOL, when should I annotate a DNA sequence vs. making sub-components? When representing a DNA sequence in SBOL as a ComponentDefinition, you can mark things like promoters and coding sequences in two different ways:

either as a SequenceAnnotation in the ComponentDefinition, or
as a Component linking to a ComponentDefinition for the sub-comoponent

When should you choose one over the other?

Typically you want to use a ComponentDefinition when an sub-structure is something you might want to pull out and re-use in another genetic design, and a SequenceAnnotation whenever it is not.

For example:

Promoters, terminators, and protein coding sequences are typically best represented with their own ComponentDefinition, since these are often pulled out and re-combined in new genetic designs.  Representing such an element with a ComponentDefinition makes it easy to find the commonalities between designs.
Assembly scars, binding sites, and cut sites are typically best represented with a SequenceAnnotation because they are typically not very interesting except in the context of their surrounding genetic context. Representing such an element with a SequenceAnnotation helps keep them associated tightly with that context.

Not every case is cut and dried, however.  For example, ribosome entry sites are often treated as separable components, but their performance is very tightly tied to the specifics of the coding sequence that they modulate.  Thus, in some cases it may make more sense to represent them separately as a ComponentDefinition and in others to represent them as a SequenceAnnotation on a joint RES/CDS ComponentDefinition.


Filtering fasta files by ID I have this code to filter some IDs from a fasta file:

for line in file_open:
    line=line.rstrip()
    if '16S' in line:
        header=line
        header=header[0] + header[1:]
        seq=file_open.readline()
        seq=seq.strip()
        print(header + '\n'  + seq)

however, this code doesn't get me all content of each ID sequence. I get this result:

>ID1
AGTTTTGCCGAATTCCTTTAGGATTTTGTGAGGAACATTGTGTTAAATGTATTTATAGTTCTACTAATTC

When I should get this:

>ID1
AGTTTTGCCGAATTCCTTTAGGATTTTGTGAGGAACATTGTGTTAAATGTATTTATAGTTCTACTAATTC
TGCGAAGGTAGCATAATCATTAGCCTTTTGATTTGAGGCTGGAATGAATGGTTTGACGAGAGATGGTCTG
GGGCCGTGGTATACTGACCATGCGAAGGTAGCATAATCATTAGCCTCTTGATTTGAGGCTGGAATGAATG
GGGCCGTGGTATAC

Thanks

You're assuming 2-line FASTA format, but the sequence can (and typically does) span multiple lines. You can use Biopython to build a simple parser:

from Bio import SeqIO

for record in SeqIO.parse("file_name.fa", "fasta"):
    if "16S" in record.description:
        print(record.format("fasta"))


IGV jump to rsID location Is it possible to load a list of rsID identifiers (e.g. rs66823261) to IGV to then be able to jump to the location of each of them when copy+pasting the rsID on the locus text box?

For human and other genomes, there is currently the feature of typing in a gene name, then jumping to the gene location.

I wonder if there is a way of achieving the same for rsIDs.

EDIT: thanks for the answers so far. I downloaded the Infinium GSA bed file from Illumina's website (https://support.illumina.com/downloads/infinium-global-screening-array-v1-0-support-files.html), loaded on an IGV session, which asked to create an index. I did a liftover of the file to hg38 and sorted with bedtools sort. Then I went to the text box, introduced one of the rs ids:

rs6569648

Got a positive result. Thank you all!

Yes, if you load the list you are interested in as a non-indexed "BED" file then you can jump to an rsID or any other name in the BED file. You'll need a targeted list, dbSNP for example is much too large to load in this way.


Genotyping with Illumina HumanOmni1-Quad before whole genome sequencing I read a paper about  (whole genome sequencing analysis of 100 Southeast Asian Malays. It says "Prior to library preparation, each sample was genotyped on the Illumina HumanOmni1-Quad as an initial DNA QC and for in-house assessment of gender and familial relations". I can hardly understand how Illumina HumanOmni1-Quad worked here for initial DNA QC and for in-house assessment of gender and familial relations. Could you please kindly explain to me a bit about this sentence?

Many thank in advance

The Illumina HumanOmni1-Quad beadchip is a microarray device consisting of 1,140,419 markers which have derived from the 1,000 genomes project. The markers chosen are is high-value regions of the genome: ABO blood typing SNPs, cSNPs, disease-associated SNPs, eSNPs, SNPs in mRNA splice sites, ADME genes, AIMs, HLA complexes, indels, introns, MHC regions, miRNA binding sites, mitochondrial DNA, PAR, promoter regions, and Y-chromosome.

The markers come together to form a sort of individual fingerprint. Germline SNPs are heritable and so can be traced up the line. Some SNPs are more prominent in subpopulations, something that Ancestry.com and other similar services leverage to predict ethnicity.

This paper utilized genomic DNA extracted from the Singapore BioBank. BioBanks tend to record age, lifestyle, family history, and health information alongside the genomic DNA. I (conservatively) assume that the Singapore BioBank at least has familial relations and gender information tied to the DNA. My interpretation is that the ivestigators SNP genotyped using the Illumina Beadchip the genomic DNA from the bank and compared their results with the "ground truth" provided by the genetic bank.

Illumina brochure for the HumanOmni1-Quad beadchip.


Print specific columns in a matrix on the basis of sample id's in the header I have a matrix file (expression.txt) in which the first column is a gene_id and from the second column the sample id's start. This matrix has 20,000 columns with sample id's corresponding to different tissue. I want to extract the samples of one tissue. Example:  SM-5GZZ7 SM-5GIEN SM-5EGGH are samples from brain tissue. So I want to separate the list of samples [(file=brain.txt;n=5000] corresponding to brain tissue from this matrix (with bash or R)
Name    SM-5GZZ7    SM-5EGHI    SM-5EGHJ    SM-5N9CS    SM-5GIEN    SM-5EGHH    SM-5EGGH
ENSG00000223972.5   0   0   0   0   0   0   0
ENSG00000227232.5   187 109 143 251 113 139 199
ENSG00000278267.1   0   0   1   0   0   0   0
ENSG00000243485.5   1   0   0   1   0   0   0
ENSG00000237613.2   0   0   0   0   0   0   0
ENSG00000268020.3   0   1   0   0   0   1   1

The file brain.txt looks like (Number of lines in this files = 5000)
SM-5GZZ7
SM-5GIEN
SM-5EGGH
SM-5GZXC
SM-5GZXO
SM-5GZY1
SM-5EGHK
SM-5EGHL

There's a million different ways to do simple taks like this. Pick one and learn it:) Python is my go to as it can do pretty much everything you'd ever want. I saved your example as csv with no spaces called expression.txt.
import pandas as pd
df = pd.read_csv('expression.txt')
df[['SM-5GZZ7','SM-5GIEN','SM-5EGGH']]

SM-5GZZ7    SM-5GIEN    SM-5EGGH
0   0   0.0 0.0
1   187 109 143 251 113 139 199 NaN NaN
2   0   0.0 0.0
3   1   0.0 0.0
4   0   0.0 0.0
5   0   0.0 1.0


Finding out proteins that have a similar domain to a specific toxin I am trying to find out that which DNA sequences(in the below) have a similar domain to a specific toxin which is Clostridium botulinum

You can find the fragment of the DNA sequence that belong to that toxin gene below:

  gi|27867582 (fragment of the known Clostridium botuninum toxin gene)
  GTGAATCAGCACCTGGACTTTCAGATGAAAAATTAAATTTAACTATCCAAAATGATGCTT
  ATATACCAAAATATGATTCTAATGGAACAAGTGATATAGAACAACATGATGTTAATGAAC
  TTAATGTATTTTTCTATTTAGATGCACAGAAAGTGCCCGAAGGTGAAAATAATGTCAATC
  TCACCTCTTCAATTGATACAGCATTATTAGAACAACCTAAAATATATACATTTTTTTCAT
  CAGAATTTATTAATAATGTCAATAAACCTGTGCAAGCAGC

I ran a blastx job over the above sequence which translates nucleotide into protein. Here's the result.

The result tells me that. The particular DNA sequence listed above, belongs to botulinum neurotoxin type E [Clostridium botulinum]

Besides, I am given two DNA sequences:

>Sequence A
TCTATCAAGTAGATTATTAAATACTACTGCACAAAATGATTCTTACGTTCCAAAATATGA
TTCTAATGGTACAAGTGAAATAAAGGAATATACTGTTGATAAACTAAATGTATTTTTCTA
TTTATATGCACAAAAAGCTCCTGAAGGTGAAAGTGCTATAAGTTTAACTTCTTCAGTTAA
TACAGCATTATTAGATGCATCTAAAGTTTATACGTTTTTTTCTTCAGATTTTATTAATAC

> Sequence B
TCCTGGCTCAGGACGAACGCTGGCGGCGTGTGCTTAACACATGCAAGTCGAGCGATGAAG
CTTCCTTCGGGAAGTGGATTAGCGGCGGACGGGTGAGTAACACGTGGGTAACCTGCCTCA
AAGTGGGGGATAGCCTTCCGAAAGGAAGATTAATACCGCATAACATAAGAGAATCGCATG
ATTTTCTTATCAAAGATTTATTGCTTTGAGATGGACCCGCGGCGCATTAGCTAGTTGGTA

Again, I translated them into proteins. The result are listed respectively(A and B):

https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Get&RID=Z423E71C015

https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Get&RID=Z42AKVWC015

The first job tells that, Sequence A belongs to botulinum neurotoxin type F [Clostridium baratii] which means that similar domain with the above botulinum neurotoxin>

Second job tells me that SequenceB belongs to unknown protein [Clostridium botulinum NCTC 2916]. So that means it is not in the same domain with the above botulinum neurotoxin right?

A protein domain is a conserved part of a given protein sequence and (tertiary) structure that can evolve, function, and exist independently of the rest of the protein chain.

This is different from global protein similarity, so similarity to neurotoxin type F doesn't necessarily mean that sequence A has similar domain to your known toxin gene. Your blast results do include conserved domain information though.


RNA-seq: How to get new expression count after normalization I've RNA seq, Human, Paired-end data, Sample size is <40. These are aligned using STAR, RSEM processed. With RSEM I've TPM and expected counts, that is two files columns as individual IDs and row as gene names.

I'm interested to normalize gene data. With edgeR tutorial (link in the end) and few other online resources I see that after following steps there's an R object that contains norm.factors (Page 15) value for each individual.

I'm unable to wrap my head around it, now. If I'm interested to get normalized gene counts, can I go ahead and multiple each individual's norm.factor into its gene counts? For example, the expected count for IND1 for Gene1 is 100 and it's norm.factor is 0.80, can I say that the normalized gene count is 100*0.80=80?

I'm not interested to perform voomm, or differential expression analysis.

Follow up question: If I'm to sample few individuals' data from these RNA-seq should I
a) normalize all together and extract interested individuals, or,
b) extract interested individuals and then perform normalization.
I think since normalization uses per person library size, each person's normalization is independent of the others.  Can someone please correct me if I'm wrong?

https://www.bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf

If I'm interested to get normalized gene counts, can I go ahead and
  multiple each individual's norm.factor into its gene counts? For
  example, the expected count for IND1 for Gene1 is 100 and it's
  norm.factor is 0.80, can I say that the normalized gene count is
  100*0.80=80?

No; you divide, but if you poke around, you can probably find a way to get edgeR to return the normalized counts.

  I think since normalization uses per person library size, each
  person's normalization is independent of the others.

Not true for the default normalization used in edgeR, TMM.  CPM normalization is independently determined for each sample.

In general, normalize everything together, unless you think that the algorithm assumptions might be violated by including all of them together.


a package with randomization/permutation tests for genomic analyses please would you let me know, is there a package in R/python/C that has implemented randomization/permutation tests for joint genomic analyses

(an example of joint genomic analyses -- when jointly considering both GENE EXPRESSION and PROTEIN BINDING along the DNA).

The context of my question is the following :

let's consider 1000 UP_REGULATED genes with increased PROTEIN X and with a HISTONE MARK Y :

in order to show that PROTEIN X is related to HISTONE MARK Y for 1000 UP-REGULATED genes, what "controls" would you use for the comparison :

-- 1000 RANDOM GENES (and multiple randomization tests)

-- 1000 UP-REG GENES with NO PROTEIN X, NO HISTONE MARK Y

-- 1000 UP-REG GENES with PROTEIN X, and NO HISTONE MARK Y

-- 1000 UP-REG GENES with NO PROTEIN X, and with HISTONE MARK Y

-- 1000 NOT-UP-REG GENES with NO PROTEIN X, NO HISTONE MARK Y

-- 1000 NOT-UP-REG GENES with PROTEIN X, and NO HISTONE MARK Y

-- 1000 NOT-UP-REG GENES with NO PROTEIN X, and with HISTONE MARK Y

-- anything else ?

any suggestions are very welcome. thanks a lot,

bogdan

ps : 've posted a similar question on BioC website : https://support.bioconductor.org/p/128828/

I would simply look at co-variance across the entire dataset and examine values >0.8 in the first instance rather than select 1000 genes.

With your approach sampling 1000 genes can be criticised for selective sampling and the way to get around this is to perform replicates, with (bootstraps) or without replacement. You could produce e.g. 1000 replicates of 1000 genes and perform co-variance on each replicate to examine Y = X.

Python has a specific permutation/combinations module call itertools, which is different to random resampling and doesn't need replicates

from itertools import combinations
X = [1,2,3]
Y = [1,2,3]
print (list(combinations(iterable, 1000) + X + Y))

The standard way is to sample is,

from random import sample
mygenes = [1, 2, 3, 4, 5]

print(sample(mygenes,3))

Output

[2, 4, 5]

There is nuance and further funcationality ... I cannot remember is "sample" is with or without replacement (its quite a big deal)


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

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

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

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

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

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

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

  -F COST

  Align DNA queries to protein reference sequences, using the specified
  frameshift cost. A value of 15 seems to be reasonable. (As a special
  case, -F0 means DNA-versus-protein alignment without frameshifts,
  which is faster.) The output looks like this:

a score=108 s prot 2  40 + 649
FLLQAVKLQDP-STPHQIVPSP-VSDLIATHTLCPRMKYQDD s dna  8 117 + 999
FFLQ-IKLWDP\STPH*IVSSP/PSDLISAHTLCPRMKSQDN

  The \ indicates a forward shift by one nucleotide, and the / indicates
  a reverse shift by one nucleotide. The * indicates a stop codon. The
  same alignment in tabular format looks like this:

108 prot 2 40 + 649 dna 8 117 + 999 4,1:0,6,0:1,10,0:-1,19

  The "-1" indicates the reverse frameshift.

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


How can I extract gene names for a metabolic pathway from KEGG? Note: this question has also been asked on Biostars

I need the get the list of gene names involved in glycolysis (to put an example). Not manually, I need to do this in a script. Ideally with Python. I think KEGG is the proper database to do this, but any other (free) suggestion would do.

How can I do it?

I whipped this little script up using the KEGG API:

#!/usr/bin/env python3
import urllib.request
import re
import sys
pathway = 'hsa00010' # glycolysis
url = "http://rest.kegg.jp/get/" + pathway
with urllib.request.urlopen(url) as f:
    lines = f.read().decode('utf-8').splitlines()
    want = 0
    for line in lines:
        fields = line.split()
        ## The list of genes starts here
        if fields[0] == 'GENE':
            want = 1
            ## The line with GENE is different
            print(fields[2].rstrip(';'))
        ## We reached the next section of the file
        elif want == 1 and re.match('^\S', line):
            sys.exit();
        ## We're still in the list of genes
        if want == 1 and len(fields)>1:
            print(fields[1].rstrip(';'))


How to get all chemical compounds that interact with a given protein target? I want to get all Pubchem Chemical compounds (CID) that interact with a given UniProt target. I have a list of UniProt targets and I need to study the relationship between these targets and chemical compounds.

Thanks

The easiest solution to the problem of mapping chemical compounds to protein IDs is to use the ChEMBL compound ID to UniProt mapping.

Doing this with PubChem is not so straight forward in part because its downloadable data format is a hideous obfuscated and badly documented mess split over multiple compressed files which in turn contain multiple compressed files...... For a publicly funded project they should really do better.


Why are there more barcodes than GEMs in 10X chromium data? Question: Why are there more barcodes than GEMs in 10X chromium data?

Introduction

I have several 10X genomics chromium libraries made with the de novo assembly/genome protocol. The 10X chromium device is advertised to produce approximately 1,000,000 GEMs (gel-emulsion droplets). In theory each GEM should have a single barcode that is used to track the provenance of each read pair. If this is true, a single library should have approximately one million barcodes, whether or not it was overloaded or underloaded.

Results

The libraries actually have more than one million barcodes, confirmed by longranger basic. Library 1 has 2,369,949 barcodes. Library 2 has 2,311,405 barcodes. Both libraries have approximately 150,000,000 read pairs.

Why this is a problem

I'm using the 10X data for some non-standard analyses and want a high number of LinkedReads per molecule (LPM). Having multiple barcodes in a GEM will reduce the number of LPMs per barcode proportional to how many barcodes were present in the GEM.

Possible explanations

If there are more than approximately one million barcodes in a library then this leaves several possibilities:

Each GEM truly has a single barcode. Therefore the Chromium device actually produced significantly more GEMs than advertised.
Each GEM has a variable number of barcodes. Therefore the number of barcodes is independent of how many GEMs are created.
Some combination of the above two options. Some variable number of GEMs and a variable number of GEMs per barcode.
Additional barcodes are created by sequencing errors and misclassification by longranger basic. This would probably manifest in many barcodes that only have one or two linked reads total.

Deeper analysis

I don't know how to test hypotheses 1, 2, & 3 bioinformatically. Hypothesis 4 is testable though by making a histogram of read counts per barcode. To do so, we will use the output of longranger basic for Lib2.

Longranger basic outputs an interleaved fastq.gz file that contains the whitelisted barcode in the header. For example (pretending the reads are shorter than 150bp):

@J00113:232:HFWCKBBXX:4:2224:17756:14889 BX:Z:AAACACCAGACAATAC-1
AATTACCCAAAATTAATATTTGTATTAATTTTCATATTAATAAC
+
JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ
@J00113:232:HFWCKBBXX:4:2224:17756:14889 BX:Z:AAACACCAGACAATAC-1
CACAATATCAAAAACAAATCAAAAATTAAGAAAATATAACTAAG
+
AAF<7FF-JJFJJFJJJJJJAJFJJJFJJJJJFJJJJJJJJJJJ

The BX:Z:AAACACCAGACAATAC-1 tells us that the index is whitelisted (BX:Z), the index sequence (AAACACCAGACAATAC), and the library index for the analysis (-1).

So, we can count the number of occurrences of each barcode with this code:

bioawk -cfastx '{print($comment)}' barcoded.fastq.gz | sort | uniq -c > bc_counts.txt

The first column of bc_counts.txt tells us how many individual reads had that barcode. We divide this number by two to see how many read pairs had that barcode and make a histogram.

awk '{print($1)/2}' barcode_counts.txt | sort | uniq -c | sort -k2 -n | awk '{printf("%d\t%d\n", $2, $1)}' > pairs_per_bc.txt

The first column of this file has the number of read pairs in the barcode, and the second column has how many barcodes have that number of read pairs. We see that there are many barcodes that only have one or two read pairs:

1       593583
2       207192
3       79409
4       34174
5       17353
6       10662
7       7997

It looks like every barcode with fewer than 10 read pairs may be noise. How many barcodes contain at least 10 reads?

>awk '{if ($1 >= 10){sum += $2} } END{print(sum)}' pairs_per_bc.txt
1347686

Around one million.

Conclusion

This analysis doesn't give conclusive evidence for the four hypotheses above, but it seems like there is a noise-generating process in 10X data that causes there to be many barcodes with only a few readpairs.

Your initial guess is almost certainly correct. I don't know about the linked read libraries, but in the 10X single cell sequencing protocol, separating real barcodes from noise barcodes is an important and sensitive step in the analysis pipeline.

The current analysis pipeline is to look at a plot, like the one you produced and find the "knees". The plot you produced has two knees, on at around 10 and one at around 200. I'd guess that the ones below 10 are what we call "error barcodes" generated by sequencing is PCR errors from real barcodes, and the ones between 10 and 200 are what we'd call "background" barcodes - in scRNAseq they represent things like GEMs that didn't encapsulate a cell, but did encapsulate some non-cellular DNA or RNA.

As there are more than 1,000,000 possible barcodes, not all of them are used.  I believe 10X knows which sequences are used and which are not (ie the barcodes are not random), so I'm surprised that these sequences make it out of the other end of the analysis.


How do I select a proper design matrix to use in DRIMSeq? I am trying to use DRIMseq for DTU with 2 treatments on two different strains of an animal model.
How do I select a proper design matrix to use in DRIMSeq ? How I know that one is the correct one to use? And how do I see which transcripts are differentially used among conditions after the pipeline was run? I suspect it is in the plotProportions function.

Also, where do I see the logFoldChange between conditions?

Many thanks in advance as it would help me a lot in my work.

> sessionInfo()
R version 3.5.3 (2019-03-11)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

Matrix products: default

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] sva_3.30.1          BiocParallel_1.16.6 genefilter_1.64.0   mgcv_1.8-28
[5] nlme_3.1-137        edgeR_3.24.3        limma_3.38.3

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.1           compiler_3.5.3       pillar_1.4.2         bitops_1.0-6
 [5] tools_3.5.3          digest_0.6.20        zeallot_0.1.0        bit_1.1-14
 [9] annotate_1.60.1      RSQLite_2.1.1        memoise_1.1.0        tibble_2.1.3
[13] lattice_0.20-38      pkgconfig_2.0.2      rlang_0.4.0          Matrix_1.2-17
[17] DBI_1.0.0            rstudioapi_0.10      yaml_2.2.0           parallel_3.5.3
[21] IRanges_2.16.0       S4Vectors_0.20.1     vctrs_0.2.0          locfit_1.5-9.1
[25] stats4_3.5.3         bit64_0.9-7          grid_3.5.3           Biobase_2.42.0
[29] AnnotationDbi_1.44.0 survival_2.44-1.1    XML_3.98-1.20        blob_1.2.0
[33] matrixStats_0.54.0   backports_1.1.4      splines_3.5.3        BiocGenerics_0.28.0
[37] xtable_1.8-4         RCurl_1.95-4.12      crayon_1.3.4

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