text:
How to compute RPKM in R?
I have the following data of fragment counts for each gene in 16 samples:
> str(expression)
'data.frame':   42412 obs. of  16 variables:
 $ sample1 : int  4555 49 122 351 53 27 1 0 0 2513 ...
 $ sample2 : int  2991 51 55 94 49 10 55 0 0 978 ...
 $ sample3 : int  3762 28 136 321 94 12 15 0 0 2181 ...
 $ sample4 : int  4845 43 193 361 81 48 9 0 0 2883 ...
 $ sample5 : int  2920 24 104 151 50 20 32 0 0 1743 ...
 $ sample6 : int  4157 11 135 324 58 26 4 0 0 2364 ...
 $ sample7 : int  3000 19 155 242 57 12 18 2 0 1946 ...
 $ sample8 : int  5644 30 227 504 91 37 11 0 0 2988 ...
 $ sample9 : int  2808 65 247 93 272 38 1108 1 0 1430 ...
 $ sample10: int  2458 37 163 64 150 29 729 2 1 1049 ...
 $ sample11: int  2064 30 123 51 142 23 637 0 0 1169 ...
 $ sample12: int  1945 63 209 40 171 41 688 3 2 749 ...
 $ sample13: int  2015 57 432 82 104 47 948 4 0 1171 ...
 $ sample14: int  2550 54 177 59 201 36 730 0 0 1474 ...
 $ sample15: int  2425 90 279 73 358 34 1052 3 3 1027 ...
 $ sample16: int  2343 56 365 67 161 43 877 3 1 1333 ...

How do I compute RPKM values from these?

What did you tried to solve this question? Did you visit Bioconductor?

First off,
Don’t use RPKMs.
They are truly deprecated because they’re confusing once it comes to paired-end reads. If anything, use FPKMs, which are mathematically the same but use a more correct name (do we count paired reads separately? No, we count fragments).
Even better, use TPM (= transcripts per million), or an appropriate cross-library normalisation method. TMP is defined as:
$$
\text{TPM}_\color{orchid}i =
    {\color{dodgerblue}{\frac{x_\color{orchid}i}{{l_\text{eff}}_\color{orchid}i}}}
    \cdot
    \frac{1}{\sum_\color{tomato}j \color{dodgerblue}{\frac{x_\color{tomato}j}{{l_\text{eff}}_\color{tomato}j}}}
    \cdot
    \color{darkcyan}{10^6}
$$
where

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

That said, FPKM an be calculated in R as follows. Note that most of the calculation happens in log transformed number space, to avoid numerical instability:
fpkm = function (counts, effective_lengths) {
    exp(log(counts) - log(effective_lengths) - log(sum(counts)) + log(1E9))
}

Here, the effective length is the transcript length minus the mean fragment length plus 1; that is, all the possible positions of an average fragment inside the transcript, which equals the number of all distinct fragments that can be sampled from a transcript.
This function handles one library at a time. I (and others) argue that this is the way functions should be written. If you want to apply the code to multiple libraries, nothing is easier using ‹dplyr›:
tidy_expression = tidy_expression %>%
    group_by(Sample) %>%
    mutate(FPKM = fpkm(Count, col_data$Lengths))

However, the data in the question isn’t in tidy data format, so we first need to transform it accordingly using ‹tidyr›:
tidy_expression = expression %>%
    gather(Sample, Count)

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

And I mentioned that TPMs are superior, so here’s their function as well:
tpm = function (counts, effective_lengths) {
    rate = log(counts) - log(effective_lengths)
    exp(rate - log(sum(exp(rate))) + log(1E6))
}

Can I ask sort of a meta question?  I've started seeing '%>%' in R code recently and had never noticed it before.  What does that do exactly?

@Greg It’s a pipe. It has been around for quite a long time (though different libraries used different operator symbols; I myself used to use %|%, similar to the shell pipe) but only recently gained mainstream traction, mainly through the dplyr library.

RPKM is defined as:

RPKM =   numberOfReads / ( geneLength/1000 * totalNumReads/1,000,000 )

As you can see, you need to have gene lengths for every gene.
Let's say geneLength is a vector which have the same number of rows as your data.frame, and every value of the vector corresponds to a gene (row) in expression.
expression.rpkm <- data.frame(sapply(expression, function(column) 10^9 * column / geneLength / sum(column)))

Regarding numerical stability
It is suggested in one of the answers, that all the computations should be performed in a log-transformed scale. In my opinion there is absolutely no reason for doing that. IEEE binary64 stores a number as binary number 1.b_{51}b{50}...b_0 times 2^{e-1023}. The relative precision doesn't depend on the exponent value given that a number is in the range [~10^-308; 10^308].
In case of RPKM we can only get out of the range if total number of reads is around 10^300, which is not realistic at all.
On the bright site there is not much harm in doing computations in the log-scale either. Worst case you loose a bit of precision.

and then picking the right gene length itself can be a nightmare depending on the organism...   So maybe you should also develop what could be the strategies to pick the correct one.

@Mitra I think this is slightly out of scope. But if you have an experience in doing that, could you add an answer here? Would be really great!

@Mitra given it a second thought, this deserves a separate question. If there will be one, I'll link it from this answer.

I don't think it is out of scope: you need the gene length to compute the RPKM. Other answers bypassed this problem by using transcript or exon length which are easily retrieved from the annotation. That is not the case with gene length which you are using in your formula.

If you are planning to do a differential expression analysis, you will probably don't need the RPKM calculation.
RPK= No.of Mapped reads/ length of transcript in kb (transcript length/1000)
RPKM = RPK/total no.of reads in million (total no of reads/ 1000000)
The whole formula together:
RPKM = (10^9 * C)/(N * L)
Where,
C = Number of reads mapped to a gene
N = Total mapped reads in the experiment
L = exon length in base-pairs for a gene

If you are looking for a more visual solution (in addition to the other answers), NCI Genomic Data Commons (TCGA repository) offers a nice formula:

Where:
RCg: Number of reads mapped to the gene
RCpc: Number of reads mapped to all protein-coding genes
RCg75: The 75th percentile read count value for genes in the sample
L: Length of the gene in base pairs

As others have pointed out, FPKMs have some problems. GDC also calculates FPKM-UQ values that are upper quartile normalized. Those are recommended for cross-sample comparison and differential expression analysis.

The denominator for FPKM-UQ is still sample specific, so this will not be an appropriate normalisation for differential expression analysis (it's within sample normalisation, not between).

GDC seems to imply otherwise: https://gdc.cancer.gov/about-data/data-harmonization-and-generation/genomic-data-harmonization/high-level-data-generation/rna-seq-quantification


text:
How to calculate overlapping genes between two genome annotation versions
I have two annotations of the same genome generated with different annotation pipelines. I want to identify overlapping gene models.
An important feature of this genome is that there are many 'genes within genes', i.e. a genemodel in the intron of another genemodel. Therefore, I only want to count two genemodels as overlapping when their coding sequence exon annotations overlap.
Using something like bedtools intersect it is straightforward to calculate overlap between the gene-level annotations.
However: I am not sure how to select genes as overlapping when only their coding sequence exons (CDS features) overlap.

Why not extract the coordinates of your CDS regions from your bed/gff files, then do the bedtools intersect?

That would still only leave me with the overlapping CDSs features. In the end I want to know the genes. Why not write your comment into an answer?

Short Answer:
In my opinion, my approach would be to pull out the CDS exons and run bedtools on those.
A Few More Details:
When you pull out the exons, make sure that you assign them all IDs if the don't already have them assigned and record which IDs "belong" to which genes. Now when you get exons that overlap, you know that they are coding and you can tie them back to which genes they originate from.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/

sort: votes


text:
Single-sample vs. joint genotyping
I am trying to understand the benefits of joint genotyping and would be grateful if someone could provide an argument (ideally mathematically) that would clearly demonstrate the benefit of joint vs. single-sample genotyping.
This is what I've gathered from other resources (Biostars, GATK forums, etc.)

Joint-genotyping helps control FDR because errors from individually genotyped samples are added up, and amplified when merging call-sets (by Heng Li on https://www.biostars.org/p/10926/)

If someone understands this, can you please clarify what is the difference on the overall FDR rate between the two scenarios (again, with an example ideally)

Greater sensitivity for low-frequency variants - By sharing information across all samples, joint calling makes it possible to “rescue” genotype calls at sites where a carrier has low coverage but other samples within the call set have a confident variant at that location. (from https://software.broadinstitute.org/gatk/documentation/article.php?id=4150)

I don't understand how the presence of a confidently called variant at the same locus in another individual can affect the genotyping of an individual with low coverage. Is there some valid argument that allows one to consider reads from another person as evidence of a particular variant in a third person? What are the assumptions for such an argument? What if that person is from a different population with entirely different allele frequencies for that variant?
Having read several of the papers (or method descriptions) that describe the latest haplotype-aware SNP calling methods (HaplotypeCaller, freebayes, Platypus) the overall framework seems to be:

Establish a prior on the allele frequency distribution at a site of interest using one (or combination) of: non-informative prior, population genetics model-based prior like Wright Fisher, prior based on established variation patterns like dbSNP, ExAC, or gnomAD.

Build a list of plausible haplotypes in a region around the locus of interest using local assembly.

Select haplotype with highest likelihood based on prior and reads data and infer the locus genotype accordingly.

At which point(s) in the above procedure can information between samples be shared or pooled? Should one not trust the AFS from a large-scale resource like gnomAD much more than the distribution obtained from other samples that are nominally party of the same "cohort" but may have little to do with each other because of different ancestry, for example?
I really want to understand the justifications and benefits offered by multi-sample genotyping and would appreciate your insights.

Say you are sequencing to 2X coverage. Suppose at a site, sample S has one reference base and one alternate base. It is hard to tell if this is a sequencing error or a heterozygote. Now suppose you have 1000 other samples, all at 2X read depth. One of them has two ALT bases; 10 of them have one REF and one ALT. It is usually improbable that all these samples have the same sequencing error. Then you can assert sample S has a het. Multi-sample calling helps to increase the sensitivity of not so rare SNPs. Note that what matters here is the assumption of error independency. Ancestry only has a tiny indirect effect.
Multi-sample calling penalizes very rare SNPs, in particular singletons. When you care about variants only, this is for good. Naively combining single-sample calls yields a higher error rate. Multi-sample calling also helps variant filtering at a later stage. For example, for a sample sequenced to 30X coverage, you would not know if a site at 45X depth is caused by a potential CNV/mismapping or by statistical fluctuation. When you see 1000 30X samples at 45X depth, you can easily know you are looking at a CNV/systematic mismapping. Multiple samples enhance most statistical signals.
Older methods pool all BAMs when calling variants. This is necessary because a single low-coverage sample does not have enough data to recover hidden INDELs. However, this strategy is not that easy to massively parallelized; adding a new sample triggers re-calling, which is very expensive as well. As we are mostly doing high-coverage sequencing these days, the old problem with INDEL calling does not matter now. GATK has this new single-sample calling pipeline where you combine per-sample gVCFs at a later stage. Such sample combining strategy is perhaps the only sensible solution when you are dealing with 100k samples.
The so-called haplotype based variant calling is a separate question. This type of approach helps to call INDELs, but is not of much relevance to multi-sample calling. Also, of the three variant callers in your question, only GATK (and Scalpel which you have not mentioned) use assembly at large. Freebayes does not. Platypus does but only to a limited extent and does not work well in practice.
I guess what you really want to talk about is imputation based calling. This approach further improves sensitivity with LD. With enough samples, you can measure the LD between two positions. Suppose at position 1000, you see one REF read and no ALT reads; at position 1500, you see one REF read and two ALT reads. You would not call any SNPs at position 1000 even given multiple samples. However, when you know the two positions are strongly linked and the dominant haplotypes are REF-REF and ALT-ALT, you know the sample under investigation is likely to have a missing ALT allele. LD transfers signals across sites and enhances the power to make correct genotyping calls. Nonetheless, as we are mostly doing high-coverage sequencing nowadays, imputation based methods only have a minor effect and are rarely applied.

Thanks, a few followups (broken out over several comments):
Top paragraph sounds more like variant calling than genotyping, which is somewhat infeasible in context of large scale high-coverage studies or ongoing clinical-focused sequencing. Do you see same logic applying to genotyping i.e. conditional on there being a variant at the locus? What I'm struggling with is understanding, batch effect detection aside, why we would trust information on AFS from a somewhat arbitrary collection of samples (genetically speaking) more than a large sample resource like gnomAD?

I don't understand why " Naively combining single-sample calls yields a higher error rate. " Can you elaborate or provide an example, in the context of genotyping?

wrt. to haplotype or assembly-based calling it's simply a reflection of the fact that best methods seem to reason over a region around a locus not just the locus itself, and me trying to figure out where in such a method information from multiple samples can be used. Specifically in the context of  high coverage whole genome sequencing where multiple batches of samples are periodically arriving for analysis and there isn't necessarily a final sample size .

Top paragraph is about deciding the genotype of sample S - it's genotyping. When you know the site allele frequency (AFS is the wrong wording here) in a larger population, you have a better prior. This prior is less accurate across populations, but better than a wright fisher prior. All these theories are only useful to lowCov. For highCov, genotype likelihood has a much bigger effect than prior and cross-sample information.

For the rest, you should ask separate questions. Pooling multiple topics under this one question is hard for you, for me and for readers to comprehend.

Anyway... On assembly-based calling, I already said it has little to do with multi-sample calling. On combining single-sample calls, I already explained in the top paragraph that jointly considering multiple samples helps genotyping – which means not using the info hurts genotyping. BTW, the body of your question is mixing variant calling and genotyping (e.g. FDR is about variant calling), which is causing confusion to me and to Devon below.

The benefit to additional samples is seen in your point 1. The likelihood of making a variant call is a function of (1) the depth of coverage supporting a given variant (ignoring mapping/base quality considerations) and (2) the likelihood of that variant existing given background knowledge. With low depth and no background knowledge, poorly covered variants will be assumed to be sequencing errors. Adding more samples can just serve then to increase the background knowledge on a position.

Thanks Devon, the question specifically targets genotyping i.e. you already have to think that there is a variant there.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
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?

What is the read length you're using?  How long are the sequences?  Do you need to hit the coverage target exactly or with some probability?

I would add a few more questions to greg's. Do I understand right that you would like to simulate sequencing of the template sequence from the file? So do the sequences represent genomes? Would you like to consider amplification bias? Sequencing errors? What sequencing platform you would like to simulate?

I am not aware of any software that can do this directly, but I would split the fasta file into one sequence per file, loop over them in BASH and invoke ART the sequence simulator (or another) on each sequence.
For more information about ART, please see their paper here: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3278762/

Could you expand on this? This really isn't a very useful answer in this state. What is ART? Where can I find it? How do I use it? How does iot work? Does it introduce errors? Can we control the error rate? Also, does it really not take multifasta files? If not, looping over several thousand files in bash will be very, very slow. There may be better ways.

Also, does it really not take multifasta files? If not, looping over several thousand files in bash will be very, very slow. There may be better ways

Where did anyone state that it does not take multifasta? OP wanted a very specific coverage per fasta record.

No, I didn't mean explain it to me. I meant please edit your answer so it provides this information. We try to have answers as comprehensive as possible and written so they can stand alone without external references on the Stack Exchange sites. Comments are hard to read, easy to miss and result in the horrible mess of biostars where you have no idea where the actual answer is. Please edit your answer to clarify.

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 ;)

This python script takes a fasta file and tsv file with counts and prints the sequences in the fasta files that many times as it is specified in the tsv file (assuming the format in the question). So if bar.tsv and foo.fasta will be your files:
from Bio import SeqIO

repeat = {}
for line in open("bar.tsv"):
    seq_id, coverage = line.split()
    repeat[seq_id] = int(coverage)

for seq_record in SeqIO.parse(foo.fasta, "fasta"):
    for i in range(repeat.get(seq_record.name, 0)):
        print(">",seq_record.name,"_",i,sep='')
        print(seq_record.seq)

Could you explain what this does and how it works? It looks like you will just be printing the original sequence $coverage times, and that can't be right. The OP needs to simulate reads, not just repeat the same sequence N times.

Ouch, then I misunderstood the question. I though that the sequences are the reads that should be simulated and 40x is laterally I would like to have this sequence 40 times. I edited the answer to at least avoid confusion.

I don't know if you misunderstood, or I did. But one of us did :) Thanks for the edit, at least now the OP can see what this does and decide for themselves.

The wgsim package by Heng Li (of BWA and samtools fame) is my go-to tool for simulating Illumina reads. It doesn't provide any convenient way to simulate differential coverage across different sequences, but it shouldn't be to hard to run wgsim multiple times, generating the desired level of coverage for each sequence of interest.
I would implement a Python script to slurp up your test file, and call wgsim (using the subprocess module) for each sequence. This will probably require you to have each sequence in a separate file. :-(

Could you expand this by adding an example command showing how to take an input fasta file and produce the fastq? This feels more like (helpful) directions to where an answer can be found than an answer.

The polyester bioconductor package can do this. It says it simulates RNA-seq reads, but I don't know if that's really any different from other NGS reads.
It can use a range of error and bias models, or learn them from a dataset.

Could you expand on this a little? How would the OP use it on their files? What range of error and bias? How would we use them? What would be a reasonable default to use? Perhaps a minimal working example would be helpful.

If you're ok with some randomness you can generate reads from your sequence file using a Poisson random variable.  You'll need to do some math to figure out what value of lambda to use in order for the expected coverage at each base pair in your read to match what you set in your text file.
For example you have a sequence S of length 1,000, a read length of 50, and an insert size of 100.  For each base b in S generate a Poisson random variable p.  You will then generate p reads from base b to b+50.  Then, generate the paired read starting at b+50+100.
Again, you would have to play with it to figure out what lambda to use but this would give you basically what you want, as long as you're ok with not having exactly the coverage you're targeting for each read.

And also will be more realistic to simulate insert distance from Poisson-like distribution! The thing is this is a reinvention of a bycicle. There are a lot of simulation tools.

Wouldn't pulling from a normal distribution be better for the insertion size?  And yeah, I know there are lots of tools but this would be one way of doing it yourself if you were interested.

from my experience, normal is a bit less flexible for fitting - the density is skewed, eg http://bioinformatics.ucdavis.edu/docs/2016-june-workshop/_images/picard-insertion-size-histogram.png . So probably a Negative Binomial will be a good choice, but I never tried to fit it to the actual insert size densities...

Simulating NGS reads while controlling sequence coverage is now easy with RNFtools (from version 0.3.1). See the tutorial, especially section Sequence extraction.
Environment preparation
First, install BioConda and add the required channels. Then either install RNFtools in the default Conda environment
conda install rnftools

or create and activate a separate Conda environment (preferable)
conda create -n rnftools rnftools
source activate rnftools

Simulation
Assume that you have a reference file ref.fa and a tab-separated coverage file coverage.tsv (e.g., those from your example). Then the following RNFtools Snakefile will do the job you want:
import rnftools
import csv

rnftools.mishmash.sample("simulation_with_coverage_control", reads_in_tuple=1)

fa = "ref.fa"
tsv = "coverage.tsv"

with open(tsv) as f:
    table = csv.reader(f, delimiter='\t')
    for seqname, cov in table:

        rnftools.mishmash.DwgSim(
            fasta=fa,
            sequences=[seqname],
            coverage=float(cov),
            read_length_1=10, # quick test with supershort reads
            read_length_2=0,
        )

include: rnftools.include()
rule: input: rnftools.input()

When you save this file (Snakefile) and run snakemake, RNFtools will simulate reads using DWGsim with the coverages defined your text file, and save all the simulated reads in simulation_with_coverage_control.fq.
You can play with all the parameters. In particular, you can use a different simulator (e.g., Art-Illumina using rnftools.mishmash.ArtIllumina). See the RNFtools documentation for more information.

Another next-gen read simulation tool is gemsim.  I haven't tested it, but I would be interested if anyone has had any experience with it.

Go to ART website. Download your package, binary or source:

Uncompress the .gz file, open the folder "art_bin_MountRainier", then "examples". Open the text file "run_test_examples_illumina.sh". This file contains numerous example of samples codes. That will give you an idea on how you can use the software. Here is an example:
# 4) simulation of paired-end reads of 100bp with the mean fragment size 500 and standard deviation 10
#   using the built-in seperated quality profiles for A, C, G, and T
$art -ss HS20 -i ./testSeq.fa -o ./paired_end_sep -l 100 -f 10 -p -m 500 -s 10 -sp -sam
#convert both aln files to a bed file
../aln2bed.pl paired_end_sep.bed paired_end_sep1.aln paired_end_sep2.aln

You can split your FASTA file sequence-by-sequence using
split -a 6 -p '^>' your_file.fa seq_

and then use any existing read simulator supporting coverage (ART, DWGsim, etc.). If you want to have all the reads mixed (not ordered by the original sequence), you can use RNFtools.
Edit 1:
As @terdon pointed out, the previous command works on OS X only. An analogical one liner for Linux (but with a slightly different naming scheme using numbers rather than letters) can be
csplit -f seq_ -n 6 your_file.fa '/^>/' {*}

To make this command work also on OS X, one needs to install coreutils (e.g., using brew) and then use gcsplit instead of csplit.
Edit 2:
Once FASTA is split by sequences, the simulation becomes straightforward and many different approaches can be used. My favorite one is using GNU Parallel. Imagine that you have your coverages in a text file called covs.txt on separate lines and in the same order as the sequences in your_file.fa, e.g.,
40
30
...

Then you can simulate reads from the original sequences using DWGsim by
ls -1 seq_* | paste covs.txt - \
    | parallel -v --colsep '\t' dwgsim -1 100 -2 100 -C {1} {2} sim_{2}

and merge the obtained FASTQ files using:
cat sim_seq_*.bwa.read1.fastq > reads.1.fq
cat sim_seq_*.bwa.read2.fastq > reads.2.fq

One possible danger of this approach is that we have assumed that the number of seq_* files is the same as the number of lines in covs.txt, which might not be true (by mistake). We should check this prior to the simulation step, e.g., by:
[[ "$(ls -1 seq_* | wc -l)" == "$(cat covs.txt | wc -l)" ]] \
    || echo "Incorrent number of lines in covs.txt"

Another caveat is that the simulated reads are not in a random order (they are grouped by source sequences).

I am afraid you're not answering the question. First of all, you are using non-standard sort options. I assume you use a Mac? The -p is a BSD-only option and not available on other *nix systems AFAIK. That said, while this will probably work to split the file on a MacOS system, splitting it is trivial and there are many, many ways of doing so. The hard bit is simulating the reads and you're not actually explaining how the OP can do that.

Thanks for the edit. However, the main issue here is that your answer is not answering the question asked. You are answering "how can I split a multifasta file into many single-sequence files" but the question is about simulating reads. Splitting the file might be part of an answer, or it might not, but it certainly isn't an answer.

Oh, and for a portable way without using any of the many specialized tools for this, you can use awk. Something like awk -vRS=">" -F'\n' 'NR>1{print ">"$0 >> $1".fa"}' file.fa. Your csplit (and your split) approach is indeed simpler though.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
How to correlate two zero inflated bedgraph-like signals?
This question pertains to iCLIP, but it could just as easily be ChIP-seq or ATAC-seq or mutation frequencies.
I have iCLIP read counts across the transcriptome and I wish to know if the signals are correlated - that is, where one of them is high, is the other likely to be high.
Often when dealing with such data (e.g. iCLIP data) we know that the data is generally sparse - that is at most positions both signals are zero and this is correct, and also zero-inflated - that is many bases that "should" have a signal are missing that data. So just calculating the Spearman's correlation is likely to give an artificially low value.
What might be a way to asses the association? I should add that the aim is the assess association of binding patterns within genes, rather than (or as well as) between genes.
Things I have thought of:

Apply some sort of smoothing to the data (eg a rolling mean). Remove any bases with 0 in both samples. Compute the spearmans.
Calculate the average pairwise distance between every read in sample one and every read in sample two. Compare this to data where the reads have been randomised within genes.

In the first case removing all bases with 0 in both samples seems wrong. But if 99.99% of all bases have zero in both samples, then this seems like its necessary for Spearman.
In the second case, the result seems like it would be non-intuitive to interpret. And also calculating this would be massively computationally intensive.

Honestly I'd just use multiBigwigSummary and then plotCorrelation from deepTools for this, but I'm a bit biased. There, the idea would be to consider each gene as a unit (you could instead use bins, but I don't think that would as nicely do what you want), namely by giving the tools a BED or GTF file input. It would then calculate the average signal in each gene/transcript and you could do your spearman's correlation. Features with 0 in all samples could optionally be removed (plotCorrelation --skipZeros).
While you certainly could go the whole 9 yards and use per-base comparisons, that seems a bit overkill and I suspect that it won't really yield appreciably more information (especially when one considers the additional time overhead).

See my clarification. - I want to look at the spatial patterns within genes, not between them.

@IanSudbery Ah, good to know, my answer isn't applicable then.

Rather than working on the base level, you could probably work on say gene level counts. Kendall's tau, an ordinal association metric, can then be used as an appropriate correlation measure.
If $X$ and $Y$ are your iCLIP replicates, $i$ represents gene index and $(x_i, y_i)$ represents the number of RBP binding sites in $X$ and $Y$ respectively for the $i^{th}$ gene, Kendall's tau is defined as :
$$
\frac{\text{#(concordant pairs)} - \text{#(discordant pairs)}}{n(n-1)/2}
$$
Where any two pairs $(x_i, y_i)$ and $(x_j, y_j)$ are concordant if:

$x_i > x_j$ AND $y_i > y_j$

OR

$x_i < x_j$ AND $y_i < y_j$

Correspondingly they are discordant if:

$x_i < x_j$ AND $y_i > y_j$

OR

$x_i > x_j$ AND $y_i < y_j$

Not sure why tex rendering is off.

See my clarification: I'm interested in the patterns within genes, rather than between them. Is tau better at dealing with data where many observations are zero?

If your interest is to study the pattern within genes, does your second sample act as a replicate? I am not sure I completely understand what "within genes" would mean here.

No, I'm trying to determine if two proteins show similar within gene binding patterns.

Mathjax support hadn't been activated. It now has been so your answer renders correctly.

It depends whether you want to treat the peak intensities as binary (comparing presence/absence of peaks in the sets) or continuous (comparing the relative magnitudes of the peaks).
Binary
For starting out, a simple binary comparison may be appropriate. You can use a peak caller of your choosing to identify peaks in each sample according to your desired criteria. Then you can use a similarity metric such as the Jaccard index to quantify the level of agreement among the peaks in the two samples.
One potential obstacle is that defining the boundaries of your peaks won't be totally straightforward. For example, a peak in one sample might have 2 overlapping peaks in the other sample, one on each end. A rough solution for this is to divide the genome into bins (maybe around 100-1000 bp, depending on your desired resolution). You can treat a peak as present in a bin if more than half of the peak lies in the bin. That way, bins in one sample can be directly compared to the corresponding bins in the other sample. Obviously, this isn't the only way to do this; other appropriate methods exist too.
Continuous
If you want to treat the peak intensities as continuous, you could apply a similar binning method, taking the "score" of a bin to be the average peak intensity at positions within that bin. You could then throw away all bins with no peaks or only low-intensity peaks throughout the genome. Then you could compute the Spearman's correlation for the remaining bins. I'm guessing it will be harder to find a strong correlation for continuous intensities, because of the amount of experimental variability that is inherently present.
If, after following these steps, the Spearman's correlation is still "artificially low" as you suggested, then this is likely a problem with the underlying data, not the overall analysis; maybe your two datasets actually don't agree that well.

I was worried about artificially low values when most observations were 0. Excluding cases where observations are 0 in both signals just didn't "feel" right.

@IanSudbery I see. So you're also interested in the extent to which the "non-peaks" agree with one another. One possibility then is not to throw away the bins where both values are zero, and to report the Spearman correlation. However, since you said the peaks are quite sparse, your zero values will make it hard to see the extent of correlation between the true positives. I agree that it's not "right" to exclude the zero values in Spearman correlation scores, but if you instead use a metric for set intersections (like the Jaccard index) excluding zeroes would be okay.

Looking at the intersection solves the sparsity problem, but not the zero-inflated one: signals from similar looking tracks tend to be close to each other, but not right on top of each other. Honestly i'm starting to think your second suggest is the best there is.

That's true, but if you use the binning approach to define intersections instead of looking for strict overlap, this could help. (i.e. it counts as an intersection if the same bin contains peaks in both tracks, even if the peaks themselves don't directly overlap)

It is one of my favorite stories.
Drop a glance to StereoGene software, it for genomic track correlation, it described in a preprint.
You also can run MACS or another peak caller and estimate the correlation of two interval sets using the GenomtriCorr package.


text:
How does one construct a cladogram of intraspecies relationships?
I have SNP data from several cultivars of rice which I have used to produce alignments, but I don't think that the usual models and algorithms used for generating phylogenetic trees are appropriate, because these cultivars are not the result of speciation events and have been interbred in their histories. How can I best calculate and visualize their degree of relatedness?

You are assuming a tree topology when you want to draw a cladogram. Under this assumption, you can just use typical tree builders. Their major assumption is the tree-like relationship.

Even with inbreeding and other genetic phenomena that might mask actual evolution of these cultivars, any phylogenetic methodology would be capable of determining relationships accurately.
Try creating a Neighbour Joining tree with MEGA, which is one of the simplest methods available. This should give you enough to check the relationships of the cultivars.

If the accessions have been interbred (conflicting with a tree-like phylogeny), you could try Splitstree, which tries to estimate a tree-like phylogeny with added hybridisation events.

Wait, you are right, but trees don't work here.
This is the area of population genetics. You will generate allele frequencies for each locus. The allele frequenceies will be used to perform Hardy-Weinberg equilibrium, F statistics, Fst and Fis
H-W start by using Structure API to work out how many populations you have. This is vital information for the next two analyses:

Fis is particularly useful here because it is the inbreeding co-efficient, the nearer 1 it is the more inbreeding the data has been through.
The Fst distances can be used to generate a tree, but this is not a model of point mutations, rather a model of migration between populations. You can input the distance matrix in MEGA and make a tree from it, Fstat API (I think) will do it for you.

Pop gen is difficult to get your head around at first, but you'll get use to it.


text:
How to convert FASTA to BED
I have a FASTA file:
> Sequence_1
GCAATGCAAGGAAGTGATGGCGGAAATAGCGTTAGATGTATGTGTAGCGGTCCC...
> Sequence_2
GCAATGCAAGGAAGTGATGGCGGAAATAGCGTTAGATGTATGTGTAGCGGTCCC....
....

I want to generate a BED file for each sequence like:
Sequence_1 0 1500
Sequence_2 0 1700

The BED regions will simply be the size of the sequences.
Q: I did that before with a one-line command. I don't remember what that is, it was on Biostars. I can't find the post now. What's the simplest way to do the conversion?

Do you mean this biostars post? https://www.biostars.org/p/15476/

@Greg The post doesn't give me one-line command to do what I want. I'm 100% sure I saw it somewhere on biostars a while ago, but I don't remember the name of the program that can do it.

What about this one?  https://www.biostars.org/p/15476/#189089

@Greg It was a Python program. I don't remember the name and thus no idea how to search.

Would we need to create a new tag, maybe something like "format-conversion", "file-conversion"? As I expect there will be many more convert from "X to Y format" type of questions.

We have many excellent answers! This will be an excellent reference for future users.
I found what exactly what I was asking in my question:

https://www.biostars.org/p/191052/

$ pip install pyfaidx
$ faidx --transform bed test.fasta > test.bed

This is the one-line command I was asking. The other answers also work, but I want to accept my own answer.

Developer of pyfaidx here. I was just about to write this answer and lament that I didn't pick a more pronounceable or memorable package name.

You can do this easily with bioawk, which is a version of awk with added features facilitating bioinformatics:
bioawk -c fastx '{print $name"\t0\t"length($seq)}' test.fa

-c fastx tells the program that the data should be parsed as fasta or fastq format. This makes the $name and $seq variables available in the awk commands.

It's good practice to have your FASTA indexed, so you can leverage the .fai you are likely to already have. If not, you can just generate the index with samtools and use some awk to make your BED:
samtools faidx $fasta
awk 'BEGIN {FS="\t"}; {print $1 FS "0" FS $2}' $fasta.fai > $fasta.bed

This will maintain tab separation but you can drop the BEGIN statement to use spaces. The BED spec only requires "whitespace" for the simple BED format.

You could adapt this awk one-liner. Note that it assumes that sequence IDs are not longer than 100 characters and that there is no description following the sequence ID on the header line.
cat myseqs.fasta | awk '$0 ~ ">" {print c; c=0;printf substr($0,2,100) "\t0\t"; } $0 !~ ">" {c+=length($0);} END { print c; }'

Otherwise, any Bio* library (Perl, Python, Ruby) provides FASTA format parsers which will extract sequence IDs and lengths.
I'd point out that whilst this resembles BED it is not, strictly-speaking, since BED maps to coordinates on a chromosome or some longer sequence object.

Inspired by this answer to a related question on read length distributions, you could do this with Biopython:
from Bio.SeqIO import parse
with open("regions.bed", "w") as bed:
    for record in parse("regions.fasta", "fasta"):
        print(record.id, 0, len(record.seq), sep="\t", file=bed)

I love Python. I prefer this answer simply because it is in Python and that is what we should be using in 2017. But why does this look so much like Perl? Readability counts. I wanted to propose a cleaner snippet, but it looks like comments do not allow code. Would it be considered good tone to edit your answer?

Point taken! I've given it a go, but you're welcome to improve on it. At present there is nobody on the site with enough reputation to accept an edit. I'm happy for you to suggest an edit through comments, or edit the post and wait until someone hits 350 rep. We can delete the meta-comments after the edit is done.

OK! I edited the answer; seems like a StackExchange-y way of doing things after all.

I too love Python, but I feel like you shouldn't need to write a script to do something as trivial as this - there are lots of command line tools that can accomplish this more quickly.

@SamStudio8 On the surface level, I agree. At the same time, if a person has a Python REPL constantly open (which is a lot of people), this has minimal overhead, and, in my humble opinion, quick algorithmic thinking should be promoted at all times in this profession. Moreover, if this conversion is part of a Snakemake pipeline, then it may be even more natural than calling shell() or writing a shell block.

@KirillG I usually have a REPL open, but I think I'd still favour using awk (or such). Years ago was a different story, but learning how to use the basics of awk has been one of the most productivity improving things I have done since switching to use Linux. But sure, I definitely agree if this is part of a larger script or pipeline, I wouldn't dare call shell. Though I'd probably use pysam over biopython, but that's just personal preference :)

@KirillG I wrote a BioPython solution that is arguably more readable/pythonic and doesn't put the whole input file into memory https://bioinformatics.stackexchange.com/a/106/104

@Chris_Rands' answer is much better than mine, please upvote it!

Thanks Scott, although with @KirillG's edit your answer is more pythonic too! However, only more recent versions of BioPython allow a string as input for SeqIO.parse (rather than a file handle) and I don't think BioPython ever explicitly closes the input file handle in this case?

Hmm. Good point about Biopython version - I have only tested the most recent version on Python 2.7.11 and 3.5.1. I was a little nervous about closing file handles too, but SeqIO.parse returns a generator, not a file, so wouldn't it automatically close the handle when we reach StopIteration? (disclaimer - I have no proof of this!)

Here is an approach with BioPython. The with statement ensures both the input and output file handles are closed and a lazy approach is taken so that only a single fasta record is held in memory at a time, rather than reading the whole file into memory, which is a bad idea for large input files. The solution makes no assumptions about the sequence ID lengths or the number of lines that the sequences are spread across:
from Bio import SeqIO

with open('sequences.fasta') as in_f, open('sequences.bed','w') as out_f:
    for record in SeqIO.parse(in_f, 'fasta'):
        out_f.write('{}\t0\t{}\n'.format(record.id, len(record)))

If the FASTA sequences are all on a single line, then the following perl one-liner should work:
cat myseqs.fasta | perl -ne 'if(/^>([^ ]+)/){print $1} else {print " 0 ",length,"\n"}'

Explanation:

If the line begins with a '>', then print everything up to the first space (but don't put a line break at the end)
Otherwise, print " 0 ", followed by the length of the line, followed by a line break


text:
Confirm success or failure of RNA-Seq normalization
I am working with a set of (bulk) RNA-Seq data collected across multiple runs, run at different times of the year. I have normalized my data using library size / quantile / RUV normalization, and would like to check (quantitatively and/or qualitatively) whether or not normalization has succeeded in removing the batch effects.
It is important to note that by "normalization has succeeded", I simply mean that unwanted variation has been removed - further analysis is required to ensure that biological variation has not been removed. What are some plots / statistical tests / software packages which provide a first-step QC for normalization?

You should use box plots and PCA plot. Let's take a look at the RUV paper:

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4404308/

Before normalization and after UQ normalization:

Libraries do not cluster as expected according to treatment. ... for UQ-normalized counts. UQ normalization does not lead to better clustering of the samples...

Before normalization, the medians in the box-plot obviously look very different among replicates.
After UQ normalization, the medians look closer but Trt.11 look like an outlier. Furthermore, the treatments aren't clustered on the PCA plot. Since they are replicates, you'd like them be close on the plot.
After RUV normalization

... RUVg shrinks the expression measures for Library 11 towards the median across libraries, suggesting robustness against outliers. ... Libraries cluster as expected by treatment. ...

The RUV has made the distribution more robust and the samples closer on the PCA plot. However, it's still not perfect as one of the treatments is not close to the other two on the first PC.
The vignettes for Bioconductor RUVSeq describes the two functions: plotRLE and plotPCA.

I might point out the distinction between an RLE plot (shown here) and an ordinary boxplot (a distinction made on page 3 of the April 2017 preprint on RLE Plots: Visualizing Unwanted Variation in High Dimensional Data at https://arxiv.org/pdf/1704.03590.pdf). Otherwise an excellent answer but important to explain what RLE means.

Visual inspection with histograms, boxplots, or some other distribution visualization is the way to go. Prior to normalization, your abundances may look something like this.

Post-normalization, they should look something like this.

See this blog post for example code.

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Improve a reference genome with sequencing data
I have a DNA sample which I know doesn't quite match my reference genome - my culture comes from a subpopulation which has undergone significant mutation since the reference was created.
The example I have in mind is E.coli. We've tried assembly using a couple of different tools and the de-novo assembly isn't as high quality as we would like, despite having tonnes of data. Approaching this from a Bayesian point of view, the reference genome provides a very good prior if we could use it wisely
From visual inspection with IGV, a significant number of both SNPs and SVs appear to be present, but an assembly built entirely from my own sequencing data is not high enough quality for my purposes.
How can I modify this reference genome to match my sample with new sequencing data (preferably with Oxford Nanopore Technologies long reads, but I can also use these to scaffold short reads if necessary), taking advantage of my knowledge that the existing reference is mostly very good, without having to access the reads which were originally used to construct the reference genome?
The goal of the project isn't to determine where the SVs are, I just need a reference that accurately represents my sample in order to use the data for downstream analysis (as the training set for machine learning.) So by a high quality reference, I mean one which represents as well as possible the sample that was sequenced. To make matters worse, this may not be the one which has the highest alignment identity if there are systematic sequencing errors, as in nanopore sequencing!

Depending on the coverage of your data and the complexity of the genome, you could either reassemble the genome de novo or run a reference guided (or reference assisted) assembly. It sounds like you're leaning more towards the latter.
There are a couple of reference-guided assembly tools available: AlignGraph and Ragout. These may or may not be appropriate depending on the organism of interest and your data types. For example, these tools are very unlikely to work well on Oxford Nanopore reads which have not been error corrected using Nanopolish or Canu -correct.

One approach to this is to use whatever data you have to iteratively update the reference genome. You can keep chain files along the way so you can convert coordinates (e.g. in gff files) from the original reference to your new pseudoreference.
A simple approach might be:

Align new data to existing reference
Call variants (e.g. samtools mpileup, GATK, or whatever is best for you)
Create new reference incorporating variants from 2
Rinse and repeat (i.e. go to 1)

You can track some simple stats as you do this - e.g. the number of new variants should decrease, the number of reads mapped should increase, and the mismatch rate should decrease, with every iteration of the above loop. Once the pseudoreference stabilises, you know you can't do much more.

You can use nanopolish using the illumina reads. Also have a look at pilon.

If it is a short-read draft assembly and you have long-reads (ONT or Pacbio) run links to scaffold the genome and then run Pilon iteratively to try to polish and fill gaps using the short-reads.

Kindel (which I wrote) can infer consensus from low quality alignments of short reads to viral references, and extending it to work with single molecule reads and larger genomes is on my to-do list, though I imagine this will require some redesign.
Presumably you're dealing with a bacterial or fungal genome in this case ? I do also have a basic C++ version, but it's a long way away from being user friendly. Anyway, it may be worth a look – feel free to get in touch with any issues you encounter. I'd use Pilon as mentioned above


text:
How to quickly determine mutations in a read of a sam file?
After DNA sequencing, I generated a sam file through alignment of a fastq file.
Before using well known variant calling programs (eg. Annovar etc.), I want to pick some reads and know what kinds of mutations are there.
Is there any quick way of doing it?

You can visualize mismatches with something like IGB

For qualitative analysis, you're probably better off using something less granular like IGV or IGB. However, if you really want to look at a couple of reads:
If you're willing to ignore sequencing errors, you can inspect the CIGAR string or the MD tag, both of which give information on the alignment of a single read.
The CIGAR string gives details on insertions, deletions, clippings, matches and mismatches. From Genome Analysis Wiki,

The sequence being aligned to a reference may have additional bases that are not in the reference or may be missing bases that are in the reference. The CIGAR string is a sequence of of base lengths and the associated operation. They are used to indicate things like which bases align (either a match/mismatch) with the reference, are deleted from the reference, and are insertions that are not in the reference. For example:

RefPos:     1  2  3  4  5  6  7     8  9 10 11 12 13 14 15 16 17 18 19
Reference:  C  C  A  T  A  C  T     G  A  A  C  T  G  A  C  T  A  A  C
Read:                   A  C  T  A  G  A  A     T  G  G  C  T
With the alignment above, you get:
POS: 5
CIGAR: 3M1I3M1D5M

Most common usage of the CIGAR string uses M (match/mismatch), I (insertion), D (deletion), S (soft clipping), and H (hard clipping). Note that = (match) and X (mismatch) are available as alternatives to the less informative M, but they are less widely used.
The MD tag gives specific details on mismatches and deletions. From the SAMtools tags specification,

The MD field aims to achieve SNP/indel calling without looking at the reference. For example, a string ‘10A5^AC6’ means from the leftmost reference base in the alignment, there are 10 matches followed by an A on the reference which is different from the aligned read base; the next 5 reference bases are matches followed by a 2bp deletion from the reference; the deleted sequence is AC; the last 6 bases are matches. The MD field ought to match the CIGAR string.

Note that neither of these will give you any idea of structural variants in short reads, and neither will be particularly readable (or helpful, due to the higher error rate) in long reads.

I wrote a program, ASCIIGenome, that I find handy in cases where you want to have a quick look at genomic data. It's a genome browser for the command line.
To view only reads containing mismatches you can use the internal function awk. To filter for reads where the NM tag (number of mismatches) is >0:
ASCIIGenome -fa genome.fa aln.bam
...

[h] for help: awk 'getSamTag("NM") > 0'

The view on terminal screen may look something like this:
Similarly, to get only reads containing indels you can use awk '$6 ~ "D|I"'
Hope this helps and feel free to report bugs & issues.

That's really cool. I like it.

samtools mpileup can do this quickly:
samtools mpileup -f reference.fasta -uv input.sam > variants.vcf

This will produce a VCF-formatted file containing information about what variants have been seen in the SAM file, aggregated for all the mapped reads.


text:
Bootstrap analysis - why it is called bootstrap?
Bootstrap analysis, bootstrapping etc are quite common jargons of bioinformatics and phylogenetics.
However, it is not very much clear to us, what exactly being meant by "a boot's straps".
Does it means a "comparison"? (Such as we hold the two straps of a boot (shoe) in "two hands"; may be an analogy to comparing 2 things (say left-'hand' side and right-'hand' side of an algebraic proof-task)) or it is something about sequence alignment (since we 'tie' a knot in boot (shoe) so that the 2 ends of the strap remain attached to each-other).
Some of my classmates even assumed it may be a process during computer booting (but I could not agree with that).
I've searched internet, but could not found anything very helpful. This wikipedia page (permalink) about bootstrapping, tells,

 Tall boots may have a tab, loop or handle at the top known as a bootstrap, allowing one to use fingers or a boot hook tool to help pulling the boots on. The saying "to pull oneself up by one's bootstraps" was already in use during the 19th century as an example of an impossible task. The idiom dates at least to 1834, when it appeared in the Workingman's Advocate: "It is conjectured that Mr. Murphee will now be enabled to hand himself over the Cumberland river or a barn yard fence by the straps of his boots." In 1860 it appeared in a comment on philosophy of mind: "The attempt of the mind to analyze itself an effort analogous to one who would lift himself by his own bootstraps".  Bootstrap as a metaphor, meaning to better oneself by one's own unaided efforts, was in use in 1922. This metaphor spawned additional metaphors for a series of self-sustaining processes that proceed without external help.
The term is sometimes attributed to a story in Rudolf Erich Raspe's The Surprising Adventures of Baron Munchausen, but in that story Baron Munchausen pulls himself (and his horse) out of a swamp by his hair (specifically, his pigtail), not by his bootstraps – and no explicit reference to bootstraps has been found elsewhere in the various versions of the Munchausen tales.

(In the quotation, highlighted portion seemed to tried to tell a meaning), but it wasn't also very much helpful. Could not correlate with bioinformatics. (Btw the same wikipedia article mentions computer booting).
Wikipedia article about Bootstrapping (statistics) (permalink) tells:

The bootstrap was published by Bradley Efron in "Bootstrap methods: another look at the jackknife" (1979), inspired by earlier work on the jackknife.Improved estimates of the variance were developed later. A Bayesian extension was developed in 1981. The bias-corrected and accelerated (BCa) bootstrap was developed by Efron in 1987, and the ABC procedure in 1992.

Still that is not very helpful.
However, the textbook, Introduction to Bioinformatics/ Arthur M. Lesk/ 3rd Edition/ Oxford / (Low Price edition) ; in its chapter 5 (Alignments and phylogenetic trees), section The problem of varying rates of evolution -> computational consideration (page 296); a clearer definition has been given:

Formal statistical tests, involving rerunning the calculation on subsets of the original data, are known as jackknifing and bootstrapping.

jackknifing is calculation with data sets sampled randomly from the original data. For phylogeny calculations from from multiple sequence alignments, select different subsets of the position in the alignment, and rerun the calculation. Finding that each subset gives the same phylogenetic tree lends it credibility. If each subset gives a different tree, none of them are trustworthy.
Bootstrapping is  similar to jackknifing except that the position chosen at random may include multiple copies of the same position, to form data sets of the same size as original, to preserve statistical properties of data sampling.

More clear definition of "Boot-strapping"; but it does not explain, where the relation exists with the boot (shoe) and its straps.

See the tag definition in Cross Validated "The bootstrap is a resampling method to estimate the sampling distribution of a statistic. " and the questions answered there

The booting of a computer is indeed related to bootstrapping: starting from no knowledge about the system state, it needs to work out what the system looks like, and how to set things up in an ordered fashion so that the computer becomes usable.
https://en.wikipedia.org/wiki/Booting

The method itself has nothing really to do with boots or straps, just as the jack-knife method also has nothing to do with knives. In the case of bootstrapping, the goal is to determine the accuracy of an estimate from random subsets of a population. Normally estimating something like the variance of the mean requires multiple independent samples, but bootstrapping allows you to perform estimates from a single population. This is essentially allowing the estimate to "pull itself up by its bootstraps", which is the only reasonable source for the term. See also the more general wikipedia article on bootstrapping outside of statistics.

So the answer is since "bootstrapping allows you to perform estimates from a single population", so the term like "standing on own feet" or "pull oneself up by own bootstraps" being used to indicate that. Have I now understood it correctly?

Yup, you got it exactly!

@AlwaysConfused “pull yourself up by your own bootstraps” idiomatically doesn’t mean the same as “standing on your own feet”. Rather, it’s a tongue-in-cheek reference to a paradoxical situation in the fictional tales of Baron Münchhausen, who (a famous liar) claims to have freed himself from a bog by pulling himself up by his own bootstraps. This is of course physically impossible, but “bootstrapping” (in various contexts, e.g. statistics and computing) at first glance appears to do just that (at second glance, not so much).

@KonradRudolph This clear comment enlightened very much. Thank you.

One of the tales of Baron Muenchhausen is so close that one would believe this to be the reason... but if you read them you'll see that he actually is pulling up himself (and his horse!) by pulling at his own hair, not his bootstraps.

In general, branch support values indicate how much of the total data support that branch in a proposed tree. But what do you do if you used all of the available data to infer the tree in the first place?
With Bootstrap supports, you generate new alignments (almost) ex nihilo. You calculate a trees support values, using its own data.
This is the connection to "pulling yourself up by your bootstraps": you don't "pull" on additionaly collected data (a rope or something in the metaphor), but rather on the very thing that your observation is based on (your boots).

Origin of name
Bootstrapping is a term of derision in many ways, prior widespread use in phylogenetics. The response here is also the story that I recall the name meaning 'pulling yourself up by your bootstraps'. If I remember the story of it origin correctly, it means its impossible to 'pull-yourself-up' with you bootstraps because there's nothing external to leverage.
Historic usage Thus traditionally, it was a last resort approach in statistics, or else a crude data mining approach, because a lot statistics at that time worked by comparing a result against a theoretical reference, i.e. critical probability. So for example, T-test, Chi-squared, ANOVA etc ... a theoretical null distribution is used or else created against which the observation is compared. The idea of not using any theoretical external distribution was not deemed cool in the days of critical probability. Bootstrapping does weird things to the distribution that few really understand (spent time investigating).
Method Technically bootstrapping is simply resampling with replacement. So if you've a data set of 100 nucleotides replicates of a resampled data set are constructed each of 100 nucleotides by randomly selecting 100 nucleotides (with replacement) from the original data set.  Performing tree estimate calculations over multiple replicates permits a robustness estimate of each 'node' in the tree (hypothetical common ancestor).
History in phylogenetics
With onset of phylogenetics however bootstrapping became a statistical norm for several different reasons, which are quite complex.  For example, attempting to perform critical probability across ever branch of a tree could run into statistical problems, what is a 'null distribution' in evolutionary biology (that question has since been reasonably well answered), in the days of the rows about the parsimony method - bootstrapping was well under radar. Hence bootstrapping is now mainstream, with some caveats.
I assume bootstrapping was originally supported by simulation studies, however it so established I don't even know it theoretical history.
Jack-knifing
'Jack-knifing' is named because it was deemed a rough and ready tool (a 'pen-knife') to assess a data set and a forerunner of the bootstrap. Its so old now I don't even remember the resampling strategy it uses, chopping the data set in half rings a bell - I honestly forget.
Today
These days there is no real reason to know the method, other than using it and the thresholds commonly accepted. Very little statistical development is based around bootstrapping. The underlying tenant is that each site is independent of the others and that is a widespread assumption across phylogenetics. That assumptions has been assessed (obviously) but the solutions were not implemented. Continually reassessing this assumption remains very interesting and could be productive longterm, but it doesn't negate the basic approach to treat each alignment position as an independent event to every other alignment position because its a baseline used throughout molecular phylogenetics.


text:
Large-scale gRNA design for a CRISPR screen
What are the best tools to design gRNAs in a high-throughput way for a CRISPR screen, e.g. targeting all protein-coding genes in a genome? I would like to take into account possible off-target effects, as well as to allow for flexibility in the PAM sequence, to have the possibility of using non-canonical PAMs.
I'm aware of CRISPRseek, but I'm afraid it might be quite slow to run for a genome-scale search.

I don't have enough experience with CRISPR to give a useful answer. I have one observation though. You say that you are afraid that CRISPRseek might be quite slow for a genome-scale search. Have you actually tried? It sounds like this problem should parallelize well. Presumably, you could run one gene at a time, running as many in parallel as your resources allow?

Perhaps GenomeCRISPR, I've personally never used it, but it has a well-documented API and appears you could automate it to go through the whole genome.


text:
Understanding the significance of BLAT score calculations for read fragment alignment
I'm attempting to reconstruct the read fragments that form a hit to a target, excluding/skipping gaps. BLAST's outfmt 6, and alternatives such as Diamond's M8 format only give the number of gap opens on the alignment. One can't just use the start and end positions of the alignment as indices to slice the read DNA sequence, as this will include the nucleotides that were gapped by the alignment.
PSL appeared to be the only format I could find that gave the actual start positions (and lengths) of each of the pieces that formed a hit (that is, each gap forces a new block).
The regular BLAST output is human readable and is a nightmare to parse. The XML format mentioned has a structure for each HSP, but an HSP can still include gaps:
  [...]
  <Hit_hsps>
    <Hsp>
      <Hsp_num>1</Hsp_num>
      <Hsp_bit-score>51.6</Hsp_bit-score>
      <Hsp_score>122</Hsp_score>
      <Hsp_evalue>9.8e-06</Hsp_evalue>
      <Hsp_query-from>2</Hsp_query-from>
      <Hsp_query-to>98</Hsp_query-to>
      <Hsp_hit-from>35</Hsp_hit-from>
      <Hsp_hit-to>65</Hsp_hit-to>
      <Hsp_query-frame>2</Hsp_query-frame>
      <Hsp_hit-frame>0</Hsp_hit-frame>
      <Hsp_identity>24</Hsp_identity>
      <Hsp_positive>27</Hsp_positive>
      <Hsp_gaps>1</Hsp_gaps>
      <Hsp_align-len>32</Hsp_align-len>
         <Hsp_qseq>ITAIGAGLQGPAGCEVIDAGGLLVMPGGIDTH</Hsp_qseq>
         <Hsp_hseq>IAAVGTGLE-PAGAEIIDAGGLLVMPGGIDVH</Hsp_hseq>
      <Hsp_midline>I A+G GL+ PAG E+IDAGGLLVMPGGID H</Hsp_midline>
    </Hsp>
  </Hit_hsps>
</Hit>

I'm trying to be lazy, I don't want to have to parse the Hsp_qseq and Hsp_hseq to determine exactly where my gaps are. I just want co-ordinates of where my read hits a protein. The PSL format (below) tells me this (tStarts):
match   mis-    rep.    N's     Q gap   Q gap   T gap   T gap   strand  Q               Q       Q       Q       T               T       T       T       block   blockSizes      qStarts  tStarts
        match   match           count   bases   count   bases           name            size    start   end     name            size    start   end     count
---------------------------------------------------------------------------------------------------------------------------------------------------------------
[...]
21      4       0       0       1       8       1       24      ++      <PROTEIN>   662     321     354     <READ>      101     2       101     2       8,17,   321,337,        2,50,

As an alternative, I've executed the BLAT aligner as a means to find alignments of translated DNA reads against a small protein database. I selected BLAT as it provided the easiest means to access information on the actual alignment blocks (with the smallest amount of parsing). That is, the PSL has a comma delimited field that enumerates each of the starting positions of a query to a target.
However, BLAT outputs a significantly greater number of alignments than alternative tools. This is likely due to its identity percentage threshold defaulting to 25% when emulating a blastx style query.
I can cope with this (apologies to my disks for the I/O), but for filtering purposes, the output PSL format appears to leave me nothing but raw match counts to work with.
The PSL gives me what I want, except now I don't have information on its significance. I know that the BLAT FAQ describes a Perl script for adding "web-based identity and score calculations". I can actually get identity percentages and scores from this script, but I cannot tell what the score actually represents. Is it comparable to bit-score? How can I recover something like an e-value from a PSL output file?

Just out of curiosity: can you elaborate on how BLAT 'provided the easiest means to access information on the actual alignment blocks (with the smallest amount of parsing)'? What about blast HSPs?

@holmrenser Don't you only get the HSPs by parsing the XML output? This is some laze on my part, as the HSPs in my XML described the hits from the perspective of the translated DNA, but I want to know where the alignments on the original sequence are without too much of my own error-prone interference. BLAT just gives the offset. It might be that I've overlooked something in the HSP data though...

No! Not at all! The default output format for blast is a simple text file and there are flags to get one HSP per line (no sequence, only the scores etc). You can get XML as well but that's just a way of making your life harder IMO. Please [edit] your question and clarify what the original (scientific) question you are asking is. Are you looking for exact matches? Distant homologs? Finally, it would be better if you could include an example of your PSL file if that's really what you want parsed.

@terdon Despite my dislike of XML parsing, it seems a more structured way to get out the components of a hit compared to trying to infer structure from a sequence of lines. For what it's worth BLAT doesn't have flags for outputting an HSP per line.

BLAT might not, I admit I've only ever used BLAT online for simple things. BLAST, however (both NCBI and wublast), certainly does. Look at the "tabular" option for the -outfm flag in NCBI blast for example. I've spent more time than I care to remember parsing blast output, if you clarify the exact question you have, I should be able to help.

@terdon Thanks, I know what options BLAST has, but I'm asking about something else. I have updated my question to try and clarify though!

Fair enough (and sorry, didn't mean to appear pushy). I was just wondering if using another tool (like BLAST which would give you the e-values in a parseable table format) wouldn't be a simpler approach. I realize that "redo the analysis with this tool instead" is not a perfect solution though :).

@terdon Hah, no worries. A lot of bioinformatics is "damn, I need to run this with another tool". I'd happily move this to chat if you have any ideas. For what it's worth, BLAST is too slow (for the size of the data) and still doesn't give the granularity on gap information that I'm looking for :( Any thoughts?

It would appear that BLAT and exonerate (as suggested by @terdon) are the only pieces of software that provide enough information for "HSPFragments" to be parsed from their outputs, at least according to BioPython's documentation, anyway:

Most search programs only have HSPs with one HSPFragment in them,
  making these two objects inseparable. However, there are programs
  (e.g. BLAT and Exonerate) which may have more than one HSPFragment
  objects in any given HSP. If you are not using these programs, you can
  safely consider HSP and HSPFragment as a single union.

It would seem the options are sticking with BLAT and generating scores afterwards with the script described in their FAQ, or using exonerate which provides a raw score as part of its default output. What these scores actually mean I will have to investigate, but for now I think the case is closed!

I would suggest using blastx and parsing the output XML. As an example, have a look at this: http://www.perlmonks.org/?node_id=1006367.
The relevant part of the hit would be the following:
<Hit>
  <Hit_num>1</Hit_num>
  <Hit_id>gnl|BL_ORD_ID|1515029</Hit_id>
  <Hit_def>43989.cce_0262 (Cyanothece ATCC 51142)</Hit_def>
  <Hit_accession>1515029</Hit_accession>
  <Hit_len>65</Hit_len>
  <Hit_hsps>
    <Hsp>
      <Hsp_num>1</Hsp_num>
      <Hsp_bit-score>40.0466</Hsp_bit-score>
      <Hsp_score>92</Hsp_score>
      <Hsp_evalue>0.00664016</Hsp_evalue>
      <Hsp_query-from>155</Hsp_query-from>
      <Hsp_query-to>253</Hsp_query-to>
      <Hsp_hit-from>12</Hsp_hit-from>
      <Hsp_hit-to>44</Hsp_hit-to>
      <Hsp_query-frame>-1</Hsp_query-frame>
      <Hsp_hit-frame>0</Hsp_hit-frame>
      <Hsp_identity>17</Hsp_identity>
      <Hsp_positive>27</Hsp_positive>
      <Hsp_gaps>0</Hsp_gaps>
      <Hsp_align-len>33</Hsp_align-len>
      <Hsp_qseq>LRGAICSMEHIEEALGKLKDWARKLIELLLGPR</Hsp_qseq>
      <Hsp_hseq>ITGAVCLMDYLEKVLEKLRELAQKLIETLLGPQ</Hsp_hseq>
      <Hsp_midline>+ GA+C M+++E+ L KL++ A+KLIE LLGP+</Hsp_midline>
    </Hsp>
  </Hit_hsps>
</Hit>

You should be able to get the position information from the following tags: <Hsp_query-from>,<Hsp_query-to>,<Hsp_hit-from>,<Hsp_hit-to>.
To translate back to you original DNA sequence I can imagine you need the <Hsp_query-frame> tag.
And then of course you get the evalue for free: <Hsp_evalue>

An HSP is still allowed to contain gaps, though. So I'd still have to parse one of the qseq or hseq :(

Correct me if I'm wrong, but the tags I mentioned (Hsp_query-from, Hsp_query-to, Hsp_hit-from, Hsp_hit-to) are coordinates on the original sequence.

You are correct (I believe)! But the HSP can contain gaps, so I can't just take those co-ordinates as given to extract the DNA from the read between those positions (as I don't want the gapped nucleotides).

The fact that the read aligns to that position, and that the coordinates are on the original sequence means you actually can. Whether or not you want to do a gapped alignment is a different matter I would think.

Not for what I need - long story short is I'm trying to "flatten" the reads into a pseudo-consensus. I need to know exactly where (what actual position) the gaps are, so I can handle them appropriately. I can't just use the starts and ends of the read as given, because taking that subsequence will include one (or more) nucleotides that should be removed to accommodate the gap(s).


text:
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?

I'm a bit late to answer this question, but we've developed pretty much exactly the system you describe, based on Git and git-annex.
Its called DataLad, free and open source software, and currently mostly known within the neurosciences.
DataLad's core data structure are DataLad datasets: Git repositories with an optional annex for large or private data. Within a dataset, code and data can be version controlled in conjunction.
You can stick to Git and git-annex commands for handling the dataset, but DataLad also offers a core API that aims to harmonize and simplify version control operations. For example, in a dataset configured according to your preferences, a datalad save commits a file based on name (glob), file type, size, or location into either Git or git-annex (instead of git add, git annex add, and git commit). And likewise, a datalad push can push Git history and annexed data to your clones (instead of git push and git-annex sync). Everything is completely decentral, and DataLad integrates via git-annex special remote mechanism with a variety of hosting providers (S3, OSF, DropBox, ...), where you could store data if you like. A high level overview is here: handbook.datalad.org/r.html?about

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.

DataLad overcomes this issue by using Git's submodule mechanishm. Datasets can nested, and operated on recursively (to preserver a monorepo-like feel). With this mechanism, we've for example exposed the human connectome project dataset (80TB, 15 million files) to GitHub: github.com/datalad-datasets/human-connectome-project-openaccess

Wow, this sounds awesome! Is there some kind of a cleanup mechanism available? E.g., the intermediate data is not needed anymore is it possible to delete it permanently?

You could use the git-annex dropunused command for this to automatically remove data that isn't referenced anymore: https://git-annex.branchable.com/git-annex-dropunused/. Or, use datalad drop / git annex drop on any annexed files of your choice to remove them permanently.

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.

By the way you can use file hash instead of timestamp at least in biomake (doi: 10.1093/bioinformatics/btx306).

@IakovDavydov I'm aware of this, but I haven't actually tried if it works

+1 for all the modification on the raw data should be documented +1 for the file hash. I am a big fan of literate programming myself and I also took advantage of that approach to document the version of my tools/packages I used for the analysis and the hash of the raw data. Hence, when I see a result or a plot I also have its context of production.

Your question is somewhat open, but I think it could prove an interesting discussion. I don't believe in many cases it is worth storing the data you have created in git. As you've noted, it isn't designed for large files (although we have git-lfs) and it's definitely not designed for binary formats such as BAM.
I'm of the opinion that how a file was created and what has been done to it since is key. Large files that took much effort to create should be mirrored somewhere (but not necessarily in a version control system). Other less-important (or less difficult to create) files that have been lost, clobbered or otherwise tainted can be regenerated as long as you know how they came to be.
For what it's worth, I've been working on a piece of software called chitin (self described as a shell for disorganised bioinformaticians). I also wrote a long blog post on why I thought this was a necessary project for me, but the main reason was that despite my attempts to organise my filesystem and make good archives of my experiments, over time I forget what my shorthand directory names meant, or exactly what program generated which data.
chitin's goal is to automatically capture changes made to the file system during the execution of a command. It knows what commands to run to re-create a particular file, what commands have used that file and can tell you when and why that file was changed (and by who ;) too).
It's not finished (nothing ever is), but I feel that you might be going down the wrong road by wanting to store all your data and its versions when really, I think most people just want to know the commands that instigated changes. If the data history is important (and your code is well versioned), then you can simply check out any commit and execute your analysis to re-generate data.

I'm sorry but I'm downvoting this. Interesting discussion point, but not answering the OP's question.

First of all, kudos to you for taking versioning seriously. The fact that you're mindful of this issue is a good sign that you want to do responsible research!
For many bioinformatics projects, data files are so large that versioning the data directly with a tool like git is impractical. But your question is really getting at a couple of different issues.

How do I do my research reproducibly, and show full provenance for each data point and result I produce?
How do I manage and synchronize my research work across multiple machines?

The short answer:

Archive the primary data.
Place your workflow under version control.
Version checksums of large data files.
Use GitHub to synchronize your workflow between machines.

The long answer:
Archive the primary data
As far as reproducibility is concerned, what matters most is the primary data: the raw, unprocessed data you collect from the instrument. If you are analyzing data that has been published by others, then write a script that will automate the task of downloading the data from its primary official source, and place that script under version control.
If you or a lab mate or a colleague produced the data and it is not yet published, then you should already have plans for submitting it to an archive. Indeed, most journals and funding agencies now require this prior to publication. I'd even go so far as to say the data should be submitted as soon as it's collected. Scientists worry a lot about having their data stolen and their ideas scooped, but statistically speaking, getting scooped is much less likely than nobody ever touching your data or reading your paper. But if you or an advisor insists, most data archives allow you to keep data private for an extended period of time until a supporting manuscript is published.
Putting (for example) Fastq files in a git repository is a bad idea for a lot of reasons. No hosting service will support files that big, git will be very slow with files that big, but most importantly git/GitHub is not archival! Use a proper data archive!
Place your workflow under version control
Treat your raw data as read-only. Only process the raw data using scripts, and keep these scripts under version control. Vince Buffalo describes this well in his book Bioinformatics Data Skills. Check it out!
Version checksums of large data files
If there are any data files that you want to track but are too big to place under version control, compute checksums and place these under version control. Checksums are very small alphanumeric strings that are, for all practical purposes, unique for each data file. So instead of putting that 5GB trimmed Fastq file or the 7GB BAM file under version control, compute their checksums and put the checksums under version control. The checksums won't tell you the contents of your files, but they can tell you when the file contents change.
This should give full disclosure and complete provenance for every data point in your analysis. The workflow has a scripts/command for downloading the primary data, scripts/commands for processing the data, and checksums that serve as a signature to validate intermediate and final output files. With this, anyone should be able to reproduce you analysis!
Use GitHub to synchronize your workflow between machines
If your workflow is already under version control with git, it's trivial to push this to a hosting service like GitHub, GitLab, or BitBucket. Then it's just a matter of using git push and git pull to keep your code up-to-date on your various machines.

The Open Science Framework uses versioning for all files and is free to use: https://osf.io
You can integrate data or code from various sources such as github, dropbox, google drive, figshare or amazon cloud
You can also store files on their server using OSF data storage, but I do not know exactly what the file size limit is.

The way we deal with this is:

All work is done in a single filesystem mounted on the cluster
This file system is mounted on local machines via sshfs/samba (depending on the location of the current "local" machine on the network).
Code is versioned with git hub
All computation is carried out via light-weight automated pipelines. We use ruffus in combination with an in-house utility layer. The system doesn't really matter as long is it no more work to add another step to the pipeline than it would be to execute it manually.
All questionably design decisions are encoded in configuration files. These configuration files, along with a very detailed log output by the pipeline (what was run, what was the git commit of the code run, what was the time stamp of the files it was run on, etc) and the initial input files are version controlled.
The benefit of this is that code + configuration + time = output. It is not expected that the whole pipeline will be rerun everytime anything is changed, but the pipeline will tell you if something is out of date (it can use timestamps or file hashes), and it can all be run in one go before publication.
Any other analysis is carried out in juptyer notebooks. These are version controlled.

To summarise, we don't synchronise because we only ever work from one disk location even if we use multiple CPU locations. We version control:

Code
Inputs, configuration, logs
Juptyer notebooks

Log records the git commits used to produce the current outputs.

Interesting Ian, this is the sort of design I aspire to, but don't really follow into practice. I'm intrigued by this in-house layer on top of ruffus, what is it?

The trick is to make writing a pipeline task easier than a cluster job submission script. The utilitiy layer is provided by CGATPipelines (www.github.com/CGATOxford/CGATPipelines)

Using Git for version-controlling code is a good practice, but it does not lend itself well to versioning large data files.  Manually syncing data across multiple nodes is asking for trouble, you want this syncing to either be handled automatically in a managed environment, or just keep the files on a single network-attached storage device.
One tool you might want to look into is Arvados, which is designed for syncing bioinformatics data and workflows across multiple machines.  From the project website:

Arvados is a platform for storing, organizing, processing, and sharing genomic and other big data. The platform is designed to make it easier for data scientists to develop analyses, developers to create genomic web applications and IT administers to manage large-scale compute and storage genomic resources. The platform is designed to run in the cloud or on your own hardware.

I would only use Arvados if it really got much better than it was before.

@nuin: What specifically do you think was lacking and/or would make it not appropriate as a solution to OP's problem?

@woemier I haven't tried it lately, and I if I remember correctly it was a PITA to setup and run simple things. But as I said, don't know if got better.

This answer will sort of only cover the big data parts, i.e. things > 100MB ish, and only if your analysis pipeline ties in with the Python ecosystem. It will require a bit of learning
Try using quilt which is sort of like 'docker for data' (Github page).
$ pip install quilt
$ quilt install uciml/iris -x da2b6f5  #note the short hash
$ python
>>> from quilt.data.uciml import iris
# you've got data

Pros:

There doesn't seem to be a filesize limit for public packages (haven't stress tested this though)
Hashes your data to ensure reproducibility
Support for versioning and tags

Cons:

Storing private data starts from $7 per user / month up to 1TB of data
Pretty much Python only at this point, with some community support for R

More information here.


text:
probeset to probeset mappings between Affymetrix arrays
I am interested in identifying mappings between different types of Affymetrix arrays. I am aware that mappings between gene and probeset can be extracted using Ensembl's Biomart database.
Ensembe gene id ENSG00000181019 maps to
1. AFFY HG-U133_Plus_2's 210519_s_at
2. AFFY HuGene-1_0-st-v1's 8002303

Is there any way to extract mappings of probeset ids between two arrays(eg: HG-U133_Plus_2, HuGene-1_0-st-v1)?
eg: 210519_s_at and 8002303

Is your question how to find overlapping probe sets between arrays?  Or how to identify probe sets in general?  It is not clear.

I have added an example. Does this make sense now?

So basically you can retrieve two gene-id mappings from Ensembl's Biomart. Correct? What is preventing you from joining these two tables based on the Ensembl gene id? You can use merge function in R or JOIN statement in SQL.

@IakovDavydov or a Unix join would work too

The problem is: All mappings are not exactly one-to-one. Either probesets of HG-U133_Plus_2 or HuGene-1_0-st-v1 mapping to one ensembl gene id.

@Prradep if you explain what the underlying problem is that you're trying to solve, we might be able to offer better advice

You should not necessarily expect to see one-to-one mappings, see my answer below.

If your question is: can probeset IDs from different platforms be mapped to one another in a similar way as mapping probesets to genes, then the answer is: Yes. BioMart allows you to map almost anything that has an ID to anything else that has an ID.
You can use BioMart either via the web interface or programatically. A brief guide to using the web interface, for mapping HG U133 Plus 2 to HuGene 1.0:

Go to the start page
Select H. sapiens Ensembl genes for your database; Human genes (GRCh38.p10) for your dataset
Click Filters in the left-hand column
Expand "REGION", scroll down to GENE and select "Input microarray probes/probesets ID list [Max 500 advised]"
Select AFFY HG U133 PLUS 2 probe ID(s) and either copy/paste or upload a list, one per line
Click Attributes in the left-hand column
Scroll through, uncheck what you don't want to see and choose what you do, for example Gene Stable ID, AFFY HuGene 1 0 st v1 probe and AFFY HG U133 Plus 2 probe
Finally click "Results" in the menu at the top of the left-hand column

And you should see a result like this (you'll need to click the "Results" button).
You should not expect that there be a one-to-one mapping, for two reasons:

Genes have multiple transcripts and probesets are mapped to each transcript
Some transcripts have more than one probeset: in this case, the HuGene IDs 8002301 and 8002303 map to transcripts for this gene

Finally: here's a programmatic example using R/BioMart:
library(biomaRt)
mart.hs <- useMart("ensembl", "hsapiens_gene_ensembl")
results <- getBM(attributes = c("ensembl_gene_id", "affy_hugene_1_0_st_v1", "affy_hg_u133_plus_2"),
                 filters = "affy_hg_u133_plus_2",
                 values = c("210519_s_at"),
                 mart = mart.hs)

results
  ensembl_gene_id affy_hugene_1_0_st_v1 affy_hg_u133_plus_2
1 ENSG00000181019               8002301         210519_s_at
2 ENSG00000181019               8002303         210519_s_at

Instead of biomaRt, it is also possible to use the mapping databases built into Bioconductor itself, and map from probe to gene, and then from gene to probe in the second. Some R code to convert between hgu133 and hgu95 using the same probe ID provided in another:
library(hgu133plus2.db)
library(hgu95av2.db)

query_probe <- "210519_s_at"

hgu133_ensembl <- select(hgu133plus2.db, keys = query_probe, columns = "ENSEMBL")

ensembl_hgu95 <- select(hgu95av2.db, keys = hgu133_ensembl$ENSEMBL, keytype = "ENSEMBL", columns = "PROBEID")

dplyr::inner_join(hgu133_ensembl, ensembl_hgu95, by = "ENSEMBL", suffix = c(".133", ".95"))

Other answers explain why there might not be one to one mapping between the probes.
The AbsID database does conversion based on mapping the probe sequences to a genome build, and then determines mappings based on overlapping genome alignment coordinates. This is really useful if you want to be sure that two probes are actually likely measuring the same transcript.
It is dependent on both probes aligning to the genome, however.
Disclaimer: I worked in the group that provides the AbsID mapping tool.