text:
Using large databases with BLAT
I'm a computer scientist working with biologists at a small school that doesn't have dedicated bioinformatics staff. I apologize if I use incorrect terminology since I have limited bioinformatics background.
Our biology staff is studying a soil sample that has been sequenced. The resulting sequences (in the hundreds of millions) were analyzed using mg-rast, but they didn't seem to have a high confidence in the results.
I've been asked to help analyze these sequences. I setup a local instance of NCBI blast, but it was far too slow to analyze that many sequences. We moved on to using BLAT. They've requested we use the nt database from NCBI. I converted the nt FASTA file to 2bit file, splitting the files into small chunks.
gfServer seems to have limitations however. When I run gfServer, I get the following error while gfServer attempts to start.

Exceeding 4 billion bases, sorry gfServer can't handle that.

Is it possible to use the full nt database using BLAT? Is there an alternative to BLAT that can handle that many sequences?

What do they actually want to do with the sequences? See the closest match? All the matches? Something else entirely?

See the closest match for each sequence, then generate a summary statistics about the sample's environment.

The closest match can be done with STAR/hisat2/bowtie2/bwa or any other standard aligner that can handle large reference fasta files. That would prove vastly faster.

If you want to use BLAT, you could split by searches by chromosome.

Why use the entire nt database? If you're analyzing a soil sample, you already know you don't care about any mammalian sequences, for instance. You can make the search much faster if you choose a subset of nt instead of the entire thing. Especially since your biologists are probably only interested in unicellular species.

If you want a drop in placement for blast that's just much faster, try DIAMOND

I'd recommend having a go with either OneCodex (commercial / Freemium) or Centrifuge.
Chris Mason's lab has recently looked at metagenomic classification and OneCodex gave the most accurate result.
Centrifuge uses Burrows-Wheeler lookups, but hasn't yet been compared by the Mason Lab.


text:
Gene ratio as imput in limma
I have a data frame with gene-expression ratio. Is it possible to input this into limma/voom to find signinficannt gene-ratios between groups of samples?
my data:
                                   S1         S2         S3         S4         S5
hsa-let-7a-5p_hsa-let-7b-3p 80.066116 197.808511 275.322581 220.421053 242.100000
hsa-let-7a-5p_hsa-let-7b-5p  2.083889   2.393050   2.337716   1.062944   1.704930
hsa-let-7a-5p_hsa-let-7c-5p 83.517241  57.388889  85.350000  40.660194  71.205882
hsa-let-7a-5p_hsa-let-7d-3p  4.199393   5.792523   4.569058   3.193290   6.022388
hsa-let-7a-5p_hsa-let-7d-5p 32.729730  17.027473  14.490662   9.507378  22.416667

What do these ratio mean? Between what are those ratios calculated? But I'm not aware of any method that used ratios to calculate differentially expressed genes, as this is like a ratio

source: Stack Exchange

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

sort: votes


text:
How to measure or assign hydrophobicity score values to individual amino acids of a PDB structure?
I want measure the hydrophobicity of each amino acid within a PDB structure file.  Since I have the PDB file I want to consider the 3D information, rather than sequence-only measures such as GRAVY.
I am not interested in measuring the SASA (solvent accessible surface area) or RD (residue depth) - although correlated, it is not the information I require.
How can I comparatively measure the hydrophobicity of the amino acids within a PDB structure file?  Is there an algorithm, library or program to do so?

Dan is correct but it doesn't answer the question.  After much thought, I have finally figured out how to do it, only thanks to Dan's answer.
The problem is that a subsequence is one-dimensional and does not therefore include any of the additional information which is extractable from the three-dimensional structure (when available).  Obviously, having the three-dimensional information is much more informative and descriptive for building classifiers, predicting interactions, interaction stability, or other purposes.

Find all intramolecular contacts within a radius of N. for the Amino Acid subsequence of interest (such as an Interface).  I set N to, e.g. 4 Angstroms.
Calculate the descriptive statistics based on the set of hydrophobicity values for all intramolecular contact atoms or amino acids (depending on the granularity).  e.g. Make an average of all of the hydrophobicity scores from the selected scale.  (I use GRAVY, although there are many others).

To improve this further, the distance value of N can be set more intelligently than as a static value.   For example, to the same as the cavity radius.

No, this doesn't makes sense (I think).
The hydrophobicity of an AA is the hydrophobicity of the AA, in whatever context, just like it's molecular weight. It's weight is it's weight, in a structure or in isolation.
I know this seems physically unrealistic, but so is an abstract measure of 'hydrophobicity', it's just the best we can do based on partition coefficients.

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


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

What are the particularities in Drop-seq that would require a different tool/method to analyze it? (It is the first time I heard about drop-seq and I am not totally sure I understood the difference between drop-seq and sc-seq)

Compared with other scRNA-seq techniques (like Smart-seq), Drop-seq sacrifices read counts for cell counts. That's the big difference. It's also very cost-effective and has been increasing in popularity since it came out in 2015. http://www.cell.com/molecular-cell/abstract/S1097-2765(17)30049-7

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

I've seen that paper and read the abstract, but I didn't know about Dropbead. Thank for the info!

Note that there isn't any license for the package. It might not be safe to use it. (Although I am sure that if you contact the authors they will grant you permission and correct this :)

The DESCRIPTION file states GNU GPL3, but yes it would be good to have a LICENSE file

SC3 is part of the Bioconductor library for R. It is based on PCA and spectral dimensionality reductions and utilizes k-means. In the introductory paper the authors perform analysis on a Drop-seq dataset of more than 44,000 cells. SC3 outputs several useful plots, including a consensus matrix, expression panel, and a marker genes heatmap. All of these utilize the hierarchical clustering.
I can't speak from personal experience but the package seems very user friendly. Being a Bioconductor package means it has lots of support and excellent documentation. I am getting ready to use it myself for a small Drop-seq dataset of olfactory epithelial cells.

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

You can try a hierarchical clustering tool HGC (hierarchical graph clustering).
You may get cell population at different hierarchy by choosing multiple cutting heights. See the paper for more details.

The current version of Seurat and Scanpy can both handle Drop-seq data. If you want to include more low-reads cells you can loosen the filtering thresholds in either of these packages. These two packages really can handle both low-read-high-cell data and high-read-low-cell data.
This ENCODE paper used Seurat to process both C1 Fluidigm data (>4000 genes per cell) and 10X data (>1000 genes per cell only).
This human limb paper used Seurat to process Drop-seq data, although they didn't annotate the cell types well. I used Seurat to reprocess their data and it works much better than what the paper presents.
I guess the key is parameter tuning.

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

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


text:
Where can I find gene expression data on one of the cell lines in NCI-60?
I'm trying to retrieve the most expressed genes from a cell line from the NCI-60, lets say cell line NCI-H23, but I'm not sure where to find the gene expression data. My main goal is to create a network using that data and iRefIndex to be loaded into Cytoscape. Does anyone know where I can find the correct data or guide me in the right direction?

CellMiner portal has plenty of data about NCI-60 cell lines, inclusive microarray gene-expression.
Also, NCI-60 is a subset of the ~1000 cell lines covered by the GDSC project for which you have multiple types of measurements (e.g. copy-number, mutation, gene-expression).

Cosmic Cell Lines Project provides gene expression data for 1000+ cell lines (including NCI-60 panel together with NCI-H23 cell line). For each gene in a cell line, you can get Z-score and over/under-expression status calculated across all tumor samples diploid for that gene.


text:
At what stage of a transcriptome assembly is it better to perform read contaminant filter?
I'm trying to assemble a bivalve transcriptome. Since bivalves are filter feeders, their transcriptomes tend to be highly contaminated by bacteria, algae and whatnot. Since I pooled several transcriptomes, I have a high amount of reads (>2B reads). I also have a not-so-high quality reference genome.
Now, I'm not sure when to perform a BLAST-based filter on these reads:

Do I filter reads with BLAST before assembling them?
Do I assemble the transcriptome and then filter out contaminant contigs?
Since I have a reference genome, do I simply not care about contaminant reads because they will not align to it?

Having this amount of reads, I'd be inclined to first filter out those reads that don't have a match against Metazoa, and then to perform a second filter after assembly in order to retain only those contigs that match against Lophotrochozoa.
If you know resources on how to manage read contamination, feel free to share them.

First, to answer your question about mapping to a low-quality reference:
1. Mapping
For mapping, low genome contiguity (low N50) doesn't really matter. You will be using a spliced aligner and short reads, so even small portions of your reads will map if they match the reference genome. What does matter though is that you are confident that the draft genome assembly does not contain any scaffolds from contaminating species like bacteria, algae, et cetera.
2. Reasons to not BLAST your raw reads
There are a few reasons that I do not suggest starting with BLAST, but other options may better suit your needs.

BLASTing 2 billion reads will take several weeks to months even on a machine with many cores (90+).
Since you're not working with a model organism, there will be many reads that do not have a BLAST hit but still are probably derived from your bivalve. I do not suggest using BLAST to filter your reads for that reason.
Similarly, many contaminating reads will not hit to anything in NCBI databases, so you will miss many of the contaminants as well. Leaving them in for now and removing with the strategy I describe in #4, then reassembling, may work better.

3. Using multiple individuals of highly heterozygous species in assembly.
Most importantly for non-model organisms and highly heterozygous species - do not assemble a transcriptome using data from more than one individual. The number of allelic differences in highly heterozygous species will make the transcript graphs a jumbled mess and you'll end up with a bunch of partial and erroneous transcripts. You can merge transcriptomes from multiple individuals later with tools like orthoFinder.
4. A strategy to filter transcripts based on mapping and BLAST
Given #1, #2, and #3 above, I propose the following strategy to filter your reads that may prevent removing true bivalve reads and will guarantee removal of reads from known contaminating species.
I would try this:

Assemble the transcriptome of each individual separately.

This avoids misassemblies due to heterozygous sites and over-filtering by Trinity or other assemblers.

Concatenate all of the transcriptomes into one file. From this point on we will use the information from all individuals to help filter.

Use somthing like cat txome_individual1.fasta txome_individual2.fasta > all_individuals.fasta.

Map all the transcripts to the reference genome s using a spliced aligner like STAR or HISAT2.
Get a bam file of all the transcripts that didn't map to the reference genome. These may be from contaminants.

Use samtools view -f 4 -b mapped_transcripts.bam | samtools sort - > output.sorted.bam to only keep the transcripts that did not mapped to the reference. -f 4 is the magic flag here that keeps things that didn't map.

Make a list of all the transcripts that did not map to the reference...

samtools view output.sorted.bam | cut -f1 | sort | uniq > unmapped_transcripts.txt

...and get only the transcripts from your fasta file that are potential contaminants.

seqtk subseq all_transcripts.fa unmapped_transcripts.txt > unmapped.fasta

BLAST all of the unmapped.fasta and identify which ones are from contaminating species. blobtools works well for this if there are too many to keep track of.
Use seqtk again to only get the known contaminating seqs from BLAST or blobtools. Let's say you call this known_contaminants.fasta
Map all reads from all individuals to the known_contaminants.fasta.
Get all the reads that map to known contaminating transcripts, known_contaminants.fasta using:

samtools view -F 4 all_reads_to_contams.bam | cut -f1 | sort | uniq > contaminating_reads.txt

Now remove all of the contaminating reads from your individual fastq files using bbmap filterbyname.sh.
Reassemble the filtered fastq files.

This strategy has the following benefits:

Does not require weeks/months of compute time to BLAST 2 billion reads.
Ensures that all reads from the reference genome are left untouched.
Ensures that all reads from known contaminants are removed.
Leaves in reads with no hits to BLAST, which could be from unassembled parts of the reference genome.
Avoids potentially filtering out low-abundance contaminants by not inputting all 2 billion reads to the transcriptome assembler, but assembling on a per-individual basis.
Leverages contaminants from different individuals to filter reads from all individuals.

Thanks for your answer, there are some interesting ideas. I want to use BLAST mainly to filter out those reads/transcript that do have a bacterial/algal hit - on my first assembly attempt, I ended up with >20kb long transcripts that were just bacterial genomes. Do allelic differences distort the graphs even with a genome-based assembly? If I assemble single individuals first, I'll have 35 individuals - not counting the pooled larvae transcriptomes, which in themselves contain 50k individuals.

Allelic differences also distort the graph with genome-based assembly, causing many misassemblies and fragmenting the output. May I ask what sort of analysis that you're doing downstream with the transcriptomes? It will only take a small amount of time to make a shell script to assemble 35 transcriptomes compared to the amount of time it will take to discover all the incomplete and incorrect transcripts if all the reads are lumped into one assembly, IMHO.
I have updated the answer to give more specific instructions and the reasoning behind each step. Hope it helps!

Wow, that's a very thorough answer! Thank you very much, it seems a very clever strategy! I want to use these transcriptomes for a variety of analyses, GSEA, 3'-UTR motifs and other.

Hi @Dunois - I think I was being a bit dramatic above, but yes, point #3 is specifically relevant for de novo transcriptome assembly. You'll get more accurate results if a transcriptome is only assembled from reads from 1 individual. If you need different life stages or tissue types, my suggestion is to assemble the transcriptomes from different individuals independently, then remove duplicate transcripts using orthofinder. With that being said, if you have multiple libraries from different tissues from the same individual, then assemble all those together.


text:
What are the right parameters to trim a small RNA transcriptome with trimmomatic?
I'm having some problems with finding the right parameters to trim my small RNA Illumina reads (51 nt long) with Trimmomatic.
Before trimming, one of the samples (21M reads) looks like this:

So for my understanding, it's quite good and the last two issues can be easily solved with adapter trimming.
I then tried to trim it with the command
java -jar /opt/Trimmomatic-0.38/trimmomatic-0.38.jar SE -threads 24 -phred33 FemaleMito3.fastq FemaleMito3_trimmed.fastq ILLUMINACLIP:adapters.ultimate.fa:2:30:10 AVGQUAL:25 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:18

but I got the following result:

So now I have a problem with the sequence quality per tile (why?), an overall drop of the quality at the end of the reads and a warning with the sequence length distribution.
Then I tried trimming just the adapters with this command:
java -jar /opt/Trimmomatic-0.38/trimmomatic-0.38.jar SE -threads 24 -phred33 FemaleMito1.fastq FemaleMito1_noadapters.fastq ILLUMINACLIP:adapters.ultimate.fa:2:30:10 AVGQUAL:25

but again I get some new issues:

Again the per tile sequence quality issue, and now a warning about the GC content.
I didn't find anywhere some precise instruction on how to trim a small RNA
library with Trimmomatic. How should I properly trim my small RNA library with Trimmomatic?

Relevant: https://bioinformatics.stackexchange.com/a/6788/96.

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

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

Thank you for your answer. I expected FastQC not being happy with the sequence duplication level and overrepresented sequences, but I couldn't explain the other two. What about the "per tile sequence quality"? Did adapters inflated artificially the quality of those tiles? Should I do something about that (filter those reads out...) or can I leave them as is?

Yes, it appears that the adapters inflated quality a bit. Don't worry about failed clusters, there's no point in trimming that minor an issue off.


text:
What is the ICGC normalized_read_count?
I downloaded gene expression data (exp_seq) from the ICGC file browser.
For each sample and gene, the file contains a normalized_read_count.
What is that value? I couldn't find any information on the ICGC website. The values are definitly too low for TPM.

By reading this thread on seqanswers and by comparing the data to TCGA, I figured out

raw_read_count is the read count which you use as input for e.g. DESeq2. It has been estimated using RSEM
normalized_read_count is equivalent to the scaled_estimate from TCGA. This is the estimated fraction of transcripts made up by a given gene, as estimated by RSEM. Multiplying this value with 1e6 yields the TPM.

Glad you could find it yourself! Many thanks for posting the answer

Thank you for your answer. This is an essential information for a user of this data. I wonder why this information on the ICGC website itself though. I found a page describing the columns on github https://github.com/icgc-dcc/dcc-docs/blob/master/docs/dictionary/release-20/sequencing-based-gene-expression-expseq-primary-file-p.md. It does not explicitly say what normalized_read_count are. Also the page corresponds to one of the older release (20), not the current one (28).


text:
Human Cell Atlas - Data availability
A news item from 2017-10-18 on the website of the Human Cell Atlas states:

In addition, the consortium today also announced the impending release of gene expression profiles from the first one million immune cells collected under the HCA, toward an initial milestone of collecting at least 30 million cells representing several tissues and organs for the atlas’ first draft. These data, to be posted on an online repository by early November, will be freely available for researchers’ use. (https://www.humancellatlas.org/news/14)

As of now, Feb 2018, I cannot find these data anywhere.
Does anybody know where I can download the gene expression data of these 1 million immune cells or when it will be eventually published?

The preview version is now available here

There is now a post on Biostars:

We are starting to define the release process for the Human Cell
  Atlas.
The first couple of HCA releases will be primarily single-cell RNA-Seq
  and we are interested to know what you would expect our releases to
  contain
Please fill out our survey if you would like to help us define our
  releases
https://www.surveymonkey.com/r/S6Y753P

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

The first production release of the DCP is happening in June.
The first data release of HCA is targetting Q4 2018.


text:
Genome assembly from error-prone reads
I understand how to assemble genome from error-free reads. I implemented like this:

Construct directed overlap graph with reads as vertices and edges as
  maximum overlap between two vertices. Edges represent the length of
  overlapping maximum substring of one read suffix with other read
  prefix. Find the path that visits all vertices only once in overlap
  graph. With small graph we can use greedy approach : Starting from
  first pattern always choose egde of maximum weight and store shortest
  common superstring spelled by this path in resulting string. Resulting
  string spelled by path is shortest common superstring of all reads,
  that is assembled genome.

For example
AACTAG
  CTAGAT

Gives the AACTAGAT as a result
Now let's say that two error-prone reads have at most one error (mismatch):
AACTAG
  CTCGAT

As you can see second read has one substitution in third character from 'A' to 'C'.
So how assembled genome should look like for the error-prone reads? Do I take A or C to resulting string? Or I'm missing something?
Note: I used short reads for example because of simplicity. I work with much larger number of reads which are 100 nucleotides long and have at most one substitution of character. So that's a big overlap graph. I know I can use De Bruijn graph, but can I use some greedy approach (like in error-free example)?

It's common for assemblers to have a coverage threshold to make sure that there are enough reads across the same region to properly assemble a contig.
You can't work out an ideal solution from a situation where the consensus sequence is ambiguous. If your overlap were representative (e.g. a heterozygous variant in a diploid genome, based on a few tens of reads of 100bp, rather than 6bp), then the resulting assembly would have a bubble in the assembly graph at the variant point.
        A
       / \
A-A-C-T   G-A-T
       \ /
        C

A conservative assembler might produce two assembled contigs corresponding to these two paths, AACTAGAT and AACTCGAT. A greedy assembler would just pick one path and use that. Subsequent assembly steps might attempt to search for bubbles and annotate the assembled contigs with variant locations.

For all mainstream short-read overlappers, if there is a mismatch, the overlap will be missed entirely. You don't get a bubble. Usually you need at least 4 reads to see a bubble and the bubble is not single-based in string/overlap graph.

Yes, for the particular example that was given, a bubble wouldn't be formed. I was assuming that the reads were representative of something with a better coverage (e.g. overlap of 50 long sequences, with half having a variant, and half not)

For short reads, the typical and the most widely used solution is to correct away sequencing errors before assembly. You can correct errors with k-mer spectrum, a trie or multi-alignment. There are many papers on this topic.
Error correction alone won't fix all sequencing errors. Remaining errors may lead to bubbles and tips in the overlap/de Bruijn graph. You can identify them by inspecting the local topology of assembly graph.

So what would be solution from error-prone sample above?

No solution to your example. You can only identify errors with sufficient coverage.


text:
single stranded cDNA vs double-stranded cDNA
Single stranded cDNA is used by race pcr and rt-pcr while next-gen sequencing (Illumina) requires double stranded cDNA. What is the reason behind this ?

I'm voting to close this question as off-topic because the question is not about bioinformatics, but for some specific methods of the wet lab. If you have a question about how to handle these sequences by computers, I think it would be better suited for the site.

Can we migrate this to biology? It's a valid question and I think it's even somewhat relevant to bioinformatic (it's important to understand generation of the data), but it would suit better to Biology.SE I think. Anyway, could you link source of your information?

Sequencing also requires single stranded DNA, otherwise the sequencing primers can’t bind. Sure there are steps in the protocol where double stranded material exists, but that’s transient.


text:
Genebank Submission
I am not sure this kind of question is allowed here. I am asking because some of you might have already experienced this.
I submitted some virus genomes in the Genebank. The genomes are on hold and are not public yet.My PI thinks thank we should rename the isolates. I have sent two emails to the admin asking that I want to rename the isolates. I sent an Excel file containing accession number, current isolate name and new name that I want to assign. It has been almost a month and I haven't heard back. So, I don't know if this is the right approach.
What should I do now ? Should I just resubmit the genomes ? Or is there any other way to do that ?
Thanks

I usually hear back from the GenBank curators within a single day.
That said, I’d still suggest you wait: it’s unlikely that they missed both of your emails (but maybe double-check the address and whether they’re likely to be caught by a spam filter): they’re probably just swamped with other submissions.
Furthermore, it’s holiday season so they might operate at drastically decreased capacity.


text:
How to get log2 fold change of RNA-Seq data for time series experiment?
I know if there is one control and one treatment group it is pretty straight forward to interpret log 2 fold change. But, I have time course experiment. I have infected cells with viruses and I harvested cells at different time points: o hr, 6 hour, 12 hr and 24 hour. Here o hr is my non-infected control. What is the best way to show log 2 fold change for each gene in this type of time course experiment. I ultimate aim to find highly expressed genes at different time points.

If you just want the change with time then use time as a continuous variable. The log2FC will then be the change per hour (you can use the default Wald test to extract it).
For time-point specific comparisons you'll instead need to treat the time points as groups and use contrasts as appropriate.


text:
Consensus sequence from SAM or BAM file?
I am trying to perform reference-based assembly. Most of the tutorials teach how to create a bam file and view alignemnts in IGV or Tablet. But, I want a assembled genome sequence in fasta format. How can I do that ? Is there any tool or linux syntax available for this purpose ?

Cross-posted on biostars: https://www.biostars.org/p/374352/

A few years ago I used AlignGraph (paper and code) for a task like this. I gave AlignGraph my raw reads, so I'm not sure if it accepts SAM/BAM input or if you have to go back to the original Fastq reads.
I'm not sure if any related/improved tools have been published since AlignGraph. Internet searches for "reference-guided assembly" return mostly transcriptome assembly tools. It might be worth looking for articles that cite the AlignGraph paper.

It seems like AlignGraph takes reads in fasta format. Does it accept fastq format ? What tool do you use to convert fastq to fasta format ?

cat reads.fastq | paste - - - - | cut -f 1-2 | sed 's/^@/>/' | tr '\t' '\n' > reads.fasta. You can replace cat with gunzip -c if your Fastqs are compressed.


text:
What are phantom peaks in ChIP-seq?
I've seen two uses of this term that seem to refer to completely different phenomena. Is my understanding correct?

In Jain et al. "Active promoters give rise to false positive ‘Phantom Peaks’ in ChIP-seq experiments" (PMC4538825), a phantom peak is a region on the genome where lots of reads align in a control experiment where there is none of the protein being IP'd, but not in a control experiment without any antibody. The authors (speculatively) attribute this to general stickiness of transcription factors and the transcriptional machinery, which can bind with potentially any antibody and get IP'd.
The ENCODE quality metrics (link) say "a cross-correlation phantom-peak is ... observed at a strand-shift equal to the read length." They are referring to a high correlation of x[i] and y[i + r] where r is the read length, x[i] contains counts of tags aligning to the Watson strand at coordinate i, and y[i+r] contains counts of tags aligning to the Crick strand at coordinate i+r. They attribute this phenomenon to variations in mappability. (I don't understand their explanation at all, but that's a separate question.)

This is almost right, but it's missing a little context. In ChIP-seq, it's typical to shift the positive and negative reads because otherwise the peaks end up offset from each other. The amount of shift is determined by trying different shifts and measuring how well the read densities correlate between positive and negative strands. In the question, the first sense talks about a peak on the plot of read density vs genomic position. The second has different x and y axes because it's talking about the cross-correlation between strands. X is the offset between positive and negative strands, and y is the correlation observed at that offset. Look at figures 4D,E in this ENCODE paper for visuals.
Landt, S. G., Marinov, G. K., Kundaje, A., Kheradpour, P., Pauli, F., Batzoglou, S., ... & Chen, Y. (2012). ChIP-seq guidelines and practices of the ENCODE and modENCODE consortia. Genome research, 22(9), 1813-1831.
https://genome.cshlp.org/content/22/9/1813.full.pdf

i'm not sure what the best terms are for strands -- "positive"/"negative"; "watson"/"crick"?

I think positive / negative (as you have here) is reasonable, especially as it relates to gene-relative direction.

Positive/negative, sense/antisense, top/bottom, forward/reverse, plus/minus... All describe the same thing.


text:
Convert CD (cluster of differentiation) into gene symbol
I want to convert names according to the "cluster of differentiation" nomenclature into accepted gene symbols. For example, CD3 should map to "CD3a;CD3b;CD3c;CD247" and "CD45RO" should map to "PTPRC". If the surface antigen isn't coded for by a gene, it should map to "" or something similar.
I've seen Obtain CD (cluster of differentiation) marker expression information, but the Abcam table used there does not give the gene symbol -- it gives a list of alternate names. I've seen table 1 of this paper, but it contains errors -- for example, TNFRSF11 should be TNFRSF11A. I could copy this table by hand, or I could download this one, though that's annoying because it doesn't have one row per CD.
Is there a canonical (or just a convenient and reliable) way to do this?

I ended up using option 4 and parsing the columns for "Synonym" and "Previous Symbol".
https://www.genenames.org/cgi-bin/genefamilies/set/471


text:
How long does it typically take to call peaks via `phantompeakqualtools`?
I'm running the script run_spp_nodups.R from phantompeakqualtools. The input is a pair of BED 3+3 ("tagalign") files that are roughly 300M, and I'm using 8 threads. How long would you expect it to take? I am happy to provide other info to the extent that I can manage -- post me a comment!

Is there a reason you want to use phantompeakqualtools? It's not exactly commonly used.

If I were building a pipeline from scratch, I probably wouldn't. I am using it as part of Anshul Kundaje's open-source pipeline.
https://github.com/kundajelab/chipseq_pipeline
I'd rather not go in and start rewiring things on my own if it's running properly -- but it's taking long enough that I suspect that something is off. Hence the question.

Is there a specific reason that people avoid it?

It's a weird hodge-podge of cryptically written and poorly documented scripts that may or may not give reasonable output. It's basically an artifact of a time when ENCODE was producing a LOT of data and no one had any clue how to reasonably QC all of it. So the NSC/RSC metrics (the main point to phantompeakquals) were created. In hindsight, they're not terribly informative, in part due to being incredibly sensitive to settings. Just do some basic QC and call macs2. The pipeline will do just that anyway.


text:
Why do I obtain very strange results with DESeq2?
I am using DESeq2 to perform a differential expression analysis, but I obtained very strange results for some genes. For some genes, I have very high log2FoldChange with very low p-value as displayed in the following table.

When I look the normalized expression values for these genes.  I can see that these genes are expressed in only one single sample (condition1= samples from 1 to 7, condition2 = samples 8 to 13).

Could you please tell me how to avoid DESeq2 to detect these genes as differentially expressed ?
I am using to following commands:
dds                      <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, directory = directory, design= ~ condition)
dds                      <- DESeq(dds)
res                      <- results(dds, contrast=c("condition",condition1,condition2))

What are the size factors for these samples (sizeFactors(dds))? I suspect that sample10 is just really wonky and throwing everything off, so you'll probably want to just exclude it.

I have check that but the factor is in the range of other samples... I have similar situations with other samples..

Have you made any of the diagnostic plots related to cook's cutoff distance? If so, post some. If not, make some (see the here for an example).

I'm guessing here, but I suspect this is at least in part due to the fact that all your samples in condition 1 are not expressed at all.
Formally you have an infinitely large fold-change between condition1 and condition2. Most tools solve this by adding a small pseudocount to unexpressed samples, but this will still result in a very high fold change. So far for the big fold-change.
If I had to guess at the reason why this is significant I would think in terms of dispersion estimation. This Bioconductor post has some discussion with the authors of DESeq2 on how this works.
Finally: calculating differential expression on these genes in these samples will never be really informative, as they really are only expressed in one sample. The way I see it you have two options:

Exclude the sample from your analysis. Straightforward, but does not solve this problem for other samples.
Filter your genes before calculating differential expression, i.e. remove all genes that are not expressed in a large number of samples. This approach is standard in for instance sleuth. A good starting point would be to remove genes that are not expressed in 75% of your samples.


text:
How can I specify to DEseq2 to only perform comparisons on which I am interested with?
I am currently performing a large RNA-seq analysis from mice PBMCs. The dataset contains around 6,000 transcriptomic profiles and I would like to use DESeq2 to identify the sets of differentially expressed genes in the different conditions. In total, I have 100 biological stimulations, and for each stimulation I have 30 control samples and 30 samples treated with a molecule of interests (so is have: (30 controls+30 treated) * 100 stimulations = 6,000 samples).
I want to identify, for each stimulation, the sets of differentially expressed genes between control samples and treated samples. I do not want to compare samples from the different stimulations. In total, I would like to have thus 100 lists of differentially expressed genes.
I have started to use deseq2 to identify these lists but deseq2 is spending lot of time to perform comparisons on which I am not interested with (comparisons between the biological stimulations).
For now, I have a table sampleTable which looks like that:

I am using DEseq2 using the following command:
DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, directory = directory, design= ~ condition)

Could you please help me with that? How can I specify to deseq2 to not perform the comparisons between the biological stimulations, but rather between control and treated sample within each stimulation ?
Thank you and best,

I am not sure that you can do that... why are you not interested between control and treated of the different stimulations?

You can specify the exact comparisons you want in the results() function. So:
dds = DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, directory = directory, design= ~ condition)
dds = DESeq(dds)
res = results(dds, contrast("condition", "treatment1", "control1"))

Note that condition should not be, "stim001-control1, stim001-control2, etc.", but instead, "stim001control, stim001control, stim001control, stim001treatment, etc.". Don't put minus signs in your levels and don't number your samples in them (they're not in different groups.
The last command would be repeated for each of the comparisons. Note that this will still be quite slow due to the number of samples and the number of groups. With so many samples, you might just make a separate sampleTable for each of the comparisons you want to make. Each would then only contain the samples relevant for that comparison, so 60 total in each.
Alternatively, if you want to fit all of the groups at once then it's likely that limma/voom will prove to have better performance. It uses different math that happens to be quicker with large models.
BTW, if the slowness is coming from DESeqDataSetFromHTSeqCount() then you can mere everything on the command line into a single matrix. A couple lines of python or even the join command can do that quickly enough.


text:
RNA velocity: competing explanations for variable ratio of spliced to unspliced transcript
In "RNA velocity of single cells", La Manno et al. look at ratios between spliced and unspliced mRNA as a way of estimating the velocity of a cell through transcriptome space.
In checking for a universal steady-state spliced-unspliced ratio, they write "We found that the steady-state behaviour of most genes across a wide range of cell types was consistent with a single fixed [ratio between spliced and unspliced] (Extended Data Fig. 3a–c). However, 11% of genes showed distinct slopes in different subsets of tissues (Extended Data Fig. 3d, e), suggesting tissue-specific alternative splicing (Extended Data Fig. 3f) or degradation rates."
Tissue-specific alternative splicing or degradation rates certainly could account for these patterns, but how do La Manno et al. rule out that some tissues are simply not in a steady state? The figure they cite (screenshotted below) uses data from the Tabula Muris adult mouse cell encyclopedia. In an adult mouse, not all cells would not be expected to be in steady state -- for example, T-cell progenitors from the thymus and epithelial cells in many organs undergo continual replenishment.

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

Could you explain why that is that case? For example, couldn't you get two slopes as follows? Suppose that keratin 17 increases during development of multiple epithelial lineages, but is constant in non-epithelial cell types. For epithelia, there would be an excess of unspliced transcripts due to the portion of the population that is developing. The epithelia would have a line with a high slope and the rest would have a line with a low slope.

The difference between developing and differentiated is a gradient.

I think I understand. So if cells are on average less differentiated in intestinal epithelia and more differentiated in tracheal epithelia, you'd expect them not to fall a line of the same slope?

You wouldn't expect undifferentiated cells to fall on any particular line, since that's a broad spectrum of states. They may have some line that they tend to fall on, but you'd then expect to see a lot of noise between that and the "differentiated line", due to various levels of differentiated status.

In addition to Devon Ryan's data-driven answer, I've found a clue in the text. La Manno et al have not used the full Tabula Muris: the figure legend indicates that they used only terminally differentiated cells.


text:
ComBat batch correction: understanding the model
A collaborator of mine is using ComBat for some RNA-seq data. I would like to understand what it's doing, and I have a specific question about the structure of the model. Equation 2.1 of the paper reads
$$Y_{ijg} = \alpha_g + X\beta_g + \gamma_{ig} + \delta_{ig}\epsilon_{ijg}$$
where $i$ is the sample, $j$ is the batch, $g$ is a gene, $Y$ is observed expression, $\alpha$ is an intercept term, $\beta$ is a vector of coefficients, $\gamma$ is an additive batch effect, $\delta$ is a multiplicative batch effect, and $X$ is a design matrix.
Why is $X$ not indexed by $i$? If there are $n$ samples, $m$ genes, and $d$ variables in the analysis, what are the dimensions of $X$ (or $X_i$) and $\beta_g$?


text:
running a samtools command for multiple bam files from 1000 genomes project
I want to compute the depth of coverage only for specific intervals in phase 3, 1000 genomes project.I have not worked with 1000 genomes project before, so a bit unfamiliar with it. I do not want to download all of the bam files for the entire genomic region, just those intervals which I am interested in them.
I guess I can use this (I found from biostars forum)
samtools -bu view 'http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG01375/alignment/HG01375.mapped.ILLUMINA.bwa.CLM.low_coverage.20120522.bam' "2:1000000-2000000" | (...)

However, I am just a bit confused. First of all, when I try the code above directly from the command line, I get an error message "No such file or directory"! Am I doing something wrong? or should I provide a specific path?
Second, I want to run  it for the mapped files for all individuals: *.bam”. through a loop.
Can anyone help me to solve my problem?

The command you're looking for is:
samtools view -c http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG01375/alignment/HG01375.mapped.ILLUMINA.bwa.CLM.low_coverage.20120522.bam 2:1000000-2000000

For a series of intervals, you'd be better off scripting something with pysam:
import pysam
bam = pysam.AlignmentFile("some file.bam")
regions = open("your regions.bed")
for region in regions:
    chrom, start, end = region.split("\t")
     bam.count(chrom, int(start), int(end))  # Do something with the result

You might see if there are bigWig files available, they'd be faster to process.

,thanks. I could rum the sam tools command for a single file. Thanks for the pysam, I will try it. But do you know how can I run a samtools command through a loop? just to remind that I do not have bam files downloaded.

You can do for loops in bash. It can be pretty painful to script anything complicated in bash, though.

yes, I know the for loop solution. my problem is that the bam files for each individual are in separate folders. I am trying this: for file in http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data//alignment/.bam;
do /data/programs/samtools-1.3.1/samtools view -c "${file}" 2:1000000-2000000
done       but get an error message.

You can't wild-card URLs, you'll need to provide them as a list.

@Anna1364 please add all this information to your question. We can't answer you if you don't tell us all the details. Give us a list of the files you want to read and explain what you need to do with those files.


text:
How to plot a graph of Linkage disequilibrium with limited information
I have a snp information in an excel file and I used SNPstat to do some analysis. But I also need a plot of LD for these snps. But my file is not formatted as the most programs required (.ped, hapmap, Plink...). Here is a screen shot of file:
I have also chromosome numbers and rs information for this snps. Due to this data format problem I couldn't use Haploview or other programs. Can you suggest me a solution to create a LD graph?

https://academic.oup.com/bioinformatics/article/28/19/2537/288671

you need an LD matrix, which your table is not.  The LD matrix is a pairwise matrix.  See this page for an example: https://www.rdocumentation.org/packages/gaston/versions/1.4.9/topics/LD.plot


text:
GATK 4.2 using IndexFeatureFile with multiple files
I am trying to use GATK IndexFeatureFile to create ".vcf.gz.tbi" file. But I have 480 files.
For the pilot study I used:
!./gatk IndexFeatureFile -I HG00096.chrom22.bam.10061.vcf.gz -O HG00096.chrom22.bam.10061.vcf.gz.tbi

Doing this procedure for 480 files is not seem very practical. Is there any way to use all file as an input and it takes all files and gives ".tbi" files as an output?
I thought
!./gatk IndexFeatureFile -I *.vcf.gz -O *.vcf.gz.tbi

or maybe not using output argument (just input). I did not worked for me.

You could run loop over the files in your shell, running a separate command for each file. For example, in Bash (or another Bourne shell):
for file in *.vcf.gz; do ./gatk IndexFeatureFile -I ${file} -O ${file}.tbi; done

Is this approach also valid for merging vcf files? Because, normally I am using bcftools for merging. For example for 3 vcf files I use this: bcftools merge HG00096.chrom22.bam.10061.vcf.gz HG00097.chrom22.bam.10062.vcf.gz HG00101.chrom22.bam.10065.vcf.gz  -o merged.vcf

In the case of bcftools merge, I guess looping shouldn't be necessary since you can merge them with one command: bcftools merge -O z -o merged.vcf.gz HG*.vcf.gz


text:
Difference between de novo transcriptome assembly methods
I have been looking around (including read the original papers) to understand what is essentially the difference between StringTie in non-reference based mode (de novo) and Trinity de novo assembly. I understand that in the genome-guided nature of StringTie, we are supported by the genome information (i.e. taking into account the mapped reads). In addition to that, we also have an option to either run using the gene annotation file or not.
Simply put, am I correct if I say that Trinity is purely looking at the overlapping reads, then assemble them, whereas StringTie (non-reference based) is looking at the mapped reads that are either close or overlapping to each other? Is there any other thing I'm missing here (I understand it's not that simple, but I'm trying to make it as simple/intuitive as possible)?
Also if that's the case, would it be safe to say that generally non-reference based StringTie is preferred over Trinity and reference based StringTie (as gene annotation can lead to some "bias", hence preventing discovery of novel TSS, exon length or even make unnecessary duplicates with small differences)?

Your intuition is correct. stringTie is just looking at clumps of alignments and how they might relate to each other (either due to spliced alignments or proximity). Trinity is doing the more computationally difficult task of finding parts of reads that overlap with each other and trying to link them together into longer sequences.
Whether it makes sense to use something like stringTie or not depends mostly on the quality of the reference sequence. If you're working on an organism with a high-quality reference sequence (pretty much anything with a lot of papers using it) then there's no reason to bother with Trinity. If your reference genome additionally has a high-quality annotation (human, mouse, etc.) then you should use it, as anything "novel" that the program is likely to output otherwise is probably just an artefact. If you either don't have an annotation or it's known to be highly problematic then it makes sense to omit it when using stringTie.
Without seeing evidence to the contrary, I don't buy the argument that using a good annotation is actually biasing the output in a problematic way. For well-studied organisms, anything reported as "novel" these days likely doesn't exist, or is otherwise so fleeting that it's biologically irrelevant. The arguments of bias sound like something that was come up with once someone didn't get as many "novel" results as they were hoping.

Thanks for your answer. About the point we will less likely find novel things in organism with high-quality annotation, I was wondering if that's always the case as there are many recent publications highlighting novel elements, even in human. I guess you are saying that generally it's sufficient, but in the case of some specific condition or localisation, finding novel elements is still relevant?

Finding a few of them, sure, but unless the system is being severely disturbed one would only expect a handful of real novel findings (those reporting more are likely primarily false-positives at this point).

It's really a misnomer to call StringTie's non-reference based mode 'de-novo.' It's still using the reference genome sequence to guide the transcript assembly, it's just not using the reference annotation. Trinity is truly de-novo in that it it assembles the transcript from the overlap of the reads without mapping them to a reference genome sequence.
If you look at the precision numbers from the StringTie paper, you'll see that it maxes out around 80% under ideal circumstances with simulated data. This means that if you're working with a well annotated genome where most transcripts have already been discovered, the vast majority of the 'novel' transcripts you find will likely turn out to be artifacts.

Thanks for the answer. I agree for the first point. For the second point, I partially agree as I think we can always find some novel transcripts in the case of uncovered condition/localisation (as my comment to Devon's answer). I hope "vast majority" here means that there could still be a few out of thousands that are indeed biologically relevant. ;)

Happy to help. It's certainly possible that some novel transcripts will be real, especially if your experimental setup is novel, but validation is likely to be a major bottleneck. In my opinion, short reads are really only suitable for quantification of known transcripts. Transcript discovery is much more effective with long read sequencing technologies like PACBIO Iso-Seq or Oxford Nanopre.


text:
How to remove all BAM read groups from all reads (not just the header)?
I have problem with one my BAMs---it appears to have invalid read groups.
Normally when I have such a problem, I remove all the read groups from the BAM header as follows:
samtools view -H your.bam | grep -v "^@RG" | samtools reheader - original.bam > new.bam

and then I would execute AddOrReplaceReadGroups from Picard:
http://broadinstitute.github.io/picard/command-line-overview.html#AddOrReplaceReadGroups
I have tried this, and there must still be reads with a read group in the optional tags RG:Z:.
How to do I definitively remove all read groups from a BAM, including all of the reads?

You might have to manually strip those auxiliary tags off:
samtools view -h your.bam | grep -v "^@RG" | sed "s/\tRG:Z:[^\t]*//" | samtools view -bo your_fixed.bam -

The sed bit is searching for the aux tag and removing everything up to the next tab.


text:
Modifying a SAM header after removing all non-primary reads
I subset a BAM to only include primary reads using the following samtools commands:
samtools view -F 256 input.bam > input.primaryOnly.sam

Now, in order to convert this SAM into a new subsetted BAM, I would need to add a header. The way I normally do this is with Picard:
java -jar picard.jar ReplaceSamHeader \
      I=input.primaryOnly.sam \
      HEADER=input.header.sam \
      O=new.sam

whereby the header input.header.sam is from the original BAM:
samtools view -H input.bam > input.header.sam

Then I could convert the SAM into a BAM with
samtools view -S -b i new.sam > new.bam

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

might save that secondary alignments exist, when they don't.
How do I get a BAM header that makes sense given my new BAM (with only primary alignments)?

"Now, in order to convert this SAM into a new subsetted BAM, I would need to add a header." I don't understand why you do ignore the sam header at the beginning ?  why don't you just use 'samtools view -F 256 -b -o out.bam input.bam' ?

samtools flagstat looks at the statistics for the reads, not the header. The thing that would possibly change in the header would be sequence lengths (i.e. from samtools idxstats), but as long as the reference sequence is the same for your original and new BAM, that won't be a problem.

@Pierre Yes, I wasn't thinking about this clearly. Thanks

Your first command only needs a slight modification to add in -h. This will create a SAM file with a header.
samtools view -h -F 256 input.bam > headered.input.primaryOnly.sam

-h        Include the header in the output.

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

And if your end-goal is a sorted BAM file, then pipe the BAM output through samtools sort:
samtools view -b -F 256 input.bam | samtools sort > sorted.input.primaryOnly.bam

So that's fixing what I think your problem is, but what about the specific question that you have asked, i.e. how can a header be added to an existing headerless SAM file?
If you have an existing BAM file, you can use Picard (as you have mentioned), or manually piece the files together with cat, which will have a similar function (although Picard probably includes a few checks to make sure that the header is consistent with the input file):
(samtools view -H input.bam; cat input.primaryOnly.sam) > headered.input.primaryOnly.sam

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


text:
Subset smaller BAM to contain several thousand rows from multiple chromosomes
There are many cases whereby I would like to subset a BAM to create a small file in order to work with (e.g. algorithmic testing, debugging, etc.)
Normally I do the following, which will subset the BAM file.bam and keep the header
samtools view -H file.bam > header.sam
samtools view file.bam | head -n 5000 | cat header.sam - | samtools view -Sb - > file.unique.bam

In this case, I would like 5000 rows in chromosome 1 and 5000 in chromosome 2.
I could first try grepping by individual chromosome, and then combining the two SAMs
e.g. here's complete BAM with grepped chr1 and (incorrect but complete) header
samtools view -H file.bam > header.sam
samtools view file.bam | grep "chr1" | cat header.sam - | samtools view -Sb - > file.unique.bam

but then I have two problems:
(1) I may not be grepping the alignments to chromsome 2---there may exist BAM rows which contain 'chr2' but are not alignments.
(2) I think one must manually edit the header. There's probably no way around this.
Is there an easy way, Bioinformatics SO?

If you're not too hung up on exact numbers like 5000 reads then you can do that with a single samtools command:
samtools view -bo subset.bam -s 123.4 alignments.bam chr1 chr2

That will select 40% (the .4 part) of the reads (123 is a seed, which is convenient for reproducibility). The convenient part of this is that it'll keep mates paired if you have paired-end reads. For 5000 reads per chromosome just change the .4 part to a sufficiently small number.
In general you don't really need to subset the header. Some tools will perform a bit better if you do, but you'll generally get the same results regardless.

You can use SAMsift:
samsift \
   -i file.bam \
   -0 'c={"chr1":5000,"chr2":5000}' \
   -f 'c[RNAME]>0' \
   -c 'c[RNAME]-=1' \
   -m nonstop-remove

Explanation:

-i file.bam – input file
-0 'c={"chr1":5000,"chr2":5000}' – initialization (create a count-down dictionary for the chromosomes of interest)
-f 'c[RNAME]>0' – filtering criterion (is the counter for the current chromosome still >0?)
-c 'c[RNAME]-=1' – code decrementing the counter of the current chromosome (5000 -> 4999 -> ... -> 1 -> 0 -> -1 -> ...)
-m nonstop-remove – remove lines causing Python errors and don't stop (in this case, an error can be caused by accessing a non-existing counter for another chromosome, e.g., for chr3)

See the SAMsift readme for more info.

If you could add some explanation about what are those options and why they are necessary. For instance I don't understand why the argument to f is c[RNAME] and not other things

I've just added the explanation.

I'd generally recommend Devon Ryan's answer. However, if you do care about having the same number from each Chr, you could use the following python/pysam code (this will output approx 5000 from every Chr):
from pysam import AlignmentFile
from random import random

nreads = 5000

infile = AlignmentFile("mybam.bam")
outfile = AlignmentFile("outbam.bam", "wb", template=infile)

for chr in infile.reference_names:
     reads_in_chr = infile.count(chr)
     frac = float(nreads)/reads_in_chr
     count = 0
     #-- Replace this loop for the exactly 5000 first reads --#
     for read in infile.fetch(chr, multiple_iterators=True):
        if random() < frac:
            count += 1
            outfile.write(read)
     print("ouputted %i reads from %s" %(count, chr))
outfile.close()

This uses all chromosomes. If you only want to use, say chr1 and chr2 replace infile.reference_names with ["chr1","chr2"].
If you want exactly 5000 reads from each chr, but don't care if they are the first ones or not, then you could replace the inner for loop with:
for read in infile.fetch(chr, multiple_iterators=True):
    count += 1
    if count >= nreads:
        break
    outfile.write(read)

if you wanted the mates of those reads as well, you could add the following after outfile.write(read):
mate_read = infile.mate(read)
outfile.write(mate_read)

Note that this is slow.

To speed things up, infile.count(chr) can be replaced with parsing idxstats() (unless one is using CRAM files).

Thanks, I'd assumed that there would be an easy way to fetch this from the index.

"To speed things up, infile.count(chr) can be replaced with parsing idxstats() " I'm not quite following. What would the above look like?

using samjdk: http://lindenb.github.io/jvarkit/SamJdk.html
$ java -jar dist/samjdk.jar --body -e \
  'Map<String,Integer> c=new HashMap<>(); public Object apply(SAMRecord r) {int n=c.getOrDefault(r.getContig(),0);if(n>=5000) return false; c.put(r.getContig(),n+1); return r;}' \
 input.bam

It'd be nice if people down-voting answers would explain why they're doing so...

@Devon (I have seen the comment until now) I downvoted because I find that a long line of code without explaining what it does internally is usually harder to understand, the code as is would improve much (IMHO) broken in several lines and with comments. It is hard for me that I don't know java to understand what is happening

If you just want to create a truncated BAM file with a header then you can significantly simplify your original code:
(samtools view -H input.bam; samtools view input.bam | head -5000) | samtools -bo output.bam

This command avoids the intermediate file and few gratuitous intermediate command invocations, at the cost of a subshell (invoked by (…)).
As above, but with formatting:
(
    samtools view -H input.bam
    samtools view input.bam | head -5000 # (*)
) \
| samtools -bo output.bam

… this can of course be extended to filter by multiple chromosomes by replacing the line marked with (*) above by one or multiple lines that subset by chromosome name (samtools view input.bam chrx, no need for grep if you have indexed the original BAM file!).

Is it samtools -bo?


text:
What is the standard way to measure contig sequence lengths in a BAM?
What is the standard way to measure contig sequence lengths in a BAM?
My understanding is that the community would use samtools idxstats to compute this information from the corresponding index file.
Are there more precise/standard approaches?

Standard approaches would be either samtools idxstats or samtools view -H, both of which will produce the exact same results. In fact, all methods will produce the exact same result, since contig length is set by the reference you align against.

Specifically, reference length is not "computed", it is stored. In theory an aligner could store anything in the contig length header records, but generally its the length, in characters, of the FASTA entry for the reference from the genome the reads were aligned to.

@IanSudbery I see your point. Finding the length of characters in the FASTA reference is somewhat of a "computation" though, right? Would "recorded" be a better way to describe it?


text:
How to subset a VCF by chromosome and keep the header?
I would like to subset a VCF which only has chromosome 2.
The problem with using various grep commands, e.g.
grep -w '^#\|^2' my.vcf > my_new.vcf

or if there's a 'chr' prefix
grep -w '^#\|^chr2' my.vcf > my_new.vcf

is that this will remove the header.
I've been trying to use vcftools with the following command:
vcftools --vcf input.vcf  --chr 2 --out subset

which does not output a VCF subset.vcf as expected, but rather subset.log:
VCFtools - 0.1.15
(C) Adam Auton and Anthony Marcketta 2009

Parameters as interpreted:
    --vcf input.vcf
    --chr 2
    --out subset

After filtering, kept 2 out of 2 Individuals
After filtering, kept 80 out of a possible 720 Sites
Run Time = 0.00 seconds

Given that the run time is Run Time = 0.00 seconds, I guess there's an error here.
How can I subset a VCF and keep the header? Is there an option with bcftools instead maybe?

This seems to work for me with bcftools filter and the -r or -R argument.

   -r, --regions chr|chr:pos|chr:from-to|chr:from-[,...]

Comma-separated list of regions, see also -R, --regions-file. Note that -r cannot be used in combination with -R.
   -R, --regions-file FILE

Regions can be specified either on command line or in a VCF, BED, or tab-delimited file (the default). The columns of the
             tab-delimited file are:
             CHROM, POS, and, optionally, POS_TO, where positions are 1-based and inclusive. Uncompressed files are stored in memory, while
  bgzip-compressed and
             tabix-indexed region files are streamed. Note that sequence names must match exactly, "chr20" is not the same as "20". Also note
  that chromosome
             ordering in FILE will be respected, the VCF will be processed in the order in which chromosomes first appear in FILE.
  However, within chromosomes, the
             VCF will always be processed in ascending genomic coordinate order no matter what order they appear in FILE. Note that
  overlapping regions in FILE can
             result in duplicated out of order positions in the output. This option requires indexed VCF/BCF files. Note that -R cannot be
  used in combination with -r.

Here's an example:
$ bcftools filter vcf_nocomp_merge_geno98.vcf.gz -r 4 | head -n 38 | colrm 100 1000000
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##fileDate=20171222
##source=PLINKv1.90
##contig=<ID=1,length=249167691>
##contig=<ID=2,length=242695901>
##contig=<ID=3,length=197800245>
##contig=<ID=4,length=190915651>
##contig=<ID=5,length=180666277>
##contig=<ID=6,length=170877445>
##contig=<ID=7,length=159086805>
##contig=<ID=8,length=146293415>
##contig=<ID=9,length=141018424>
##contig=<ID=10,length=135434552>
##contig=<ID=11,length=134938471>
##contig=<ID=12,length=133763353>
##contig=<ID=13,length=115045730>
##contig=<ID=14,length=107285438>
##contig=<ID=15,length=102369712>
##contig=<ID=16,length=90141356>
##contig=<ID=17,length=81006630>
##contig=<ID=18,length=78014583>
##contig=<ID=19,length=59071322>
##contig=<ID=20,length=62906515>
##contig=<ID=21,length=48077813>
##contig=<ID=22,length=51156934>
##contig=<ID=23,length=154847490>
##INFO=<ID=PR,Number=0,Type=Flag,Description="Provisional reference allele, may not be based on rea
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##bcftools_filterVersion=1.2-187-g1a55e45+htslib-1.2.1-256-ga356746
##bcftools_filterCommand=filter -r 4 vcf_nocomp_merge_geno98.vcf.gz
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  1_125   3_5555  4_7565  5_9
4   71566   rs13125929  T   C   .   PASS    PR  GT  0/1 0/1 0/0
4   87113   rs2006748   T   C   .   PASS    PR  GT  0/1 0/0 0/0
4   110646  rs11727494  C   T   .   PASS    PR  GT  0/1 0/0 1/1
4   142550  rs11735742  T   C   .   PASS    PR  GT  0/1 0/0 0/0
4   200631  rs6826124   A   G   .   PASS    PR  GT  0/1 1/1 0/1
4   221623  rs7695945   G   T   .   PASS    PR  GT  0/0 0/0 0/0

You can use grep for this. The problem is the -w flag which means "match entire words", but VCF header lines start with ##, not #. Here are some alternative grep commands:
grep -E '^(#|2[[:space:]])' my.vcf

Or, if you have GNU grep (default on Linux):
grep -P '^(#|2\t)' my.vcf

To match either 2 or chr2:
grep -E '^(#|(chr2|2)[[:space:]])' my.vcf

or
grep -P '^(#|(chr){0,1}2\t)' my.vcf

Alternatively, you can use awk:
awk '/^#/ || $1=="2" || $2=="chr2"' my.vcf

Finally, if you want to use vcftools, you need to also use the --recode option:

--recode
--recode-bcf
These  options  are used to generate a new file in either VCF or BCF from the input VCF or BCF file after applying   the
  filtering options specified by the user. The output file has the
  suffix  ".recode.vcf"  or  ".recode.bcf".  By   default,  the  INFO
  fields are removed from the output file, as the INFO values may be
  invalidated by the recoding

(e.g. the total depth may need to be recalculated if individuals are removed). This behavior may be  overriden  by
  the following options. By default, BCF files are written out as BGZF compressed files.
But that will remove the INFO fields as explained above, so just use normal text-parsing tools like the ones I suggested.

grep (at least the version on Linux) can use multiple -e commands as expressions to search for; lines matching any expression will be displayed.

@gringer yes, that's POSIX. Personally, I prefer ERE or PCRE syntax over multiple -e, but both work.

These are both great answers. I will probably mark @gringer as the accepted as there are less opportunities for readers to make mistakes, but personally, I'm using these grep commands. Thank you!

@EB2127 you're very welcome, and please feel very free to mark whichever you prefer. I personally prefer rolling my own parser scripts, but with formats as complex as vcf, it is indeed usually better to use dedicated tools. I would go for gringer's answer as well in your place :)

You can do this more easily with the GATK tool SelectVariants, which allows you to subset variants based on various criteria. Here's the doc page:
https://gatk.broadinstitute.org/hc/en-us/articles/360041849831-SelectVariants
For this case (subsetting by chromosome) you would simply do:
gatk SelectVariants \
     -R reference.fasta \
     -V input.vcf \
     -L chr2 \
     -O output.chr2.vcf

If you wanted to do this for several chromosomes, or regions within chromosomes, you would give the tool a file of intervals with -L instead of the chromosome name. Suitable interval file formats are described here:
https://gatk.broadinstitute.org/hc/en-us/articles/360035531852-Intervals-and-interval-lists

In the folder of the command ( ./gatk ), the Fasta index file and the Fasta dict file of the reference.fasta should also be included.


text:
Tools for comparing/visualizing FASTAs?
I have two FASTAs/assemblies of the same species, but the bases are somewhat different. I would like to explore this.
What tools/methods exist to compare two FASTAs and see the difference in structure between them? Are there visualization tools available, something like IGV?
I think users would need to do a multiple-sequence alignment between these FASTAs, and then check for the differences between them. However, I don't know if a "standard approach" exists to do this.

You can try Mauve multiple genome aligner.

Mauve is a very good tool for this! Thanks. Does Mauve allow for us to visualize base-level differences (easily)?

I am not sure. I have not tried it at that resolution and i have not used it in a long time.

Depending on how long the sequences are, you could use the LAST Web Service:
http://lastweb.cbrc.jp/
It has an HTML view that shows base-level differences with colour-coding.
LAST [download / command-line version] can also be used for producing whole-genome comparison dot-plots.
Another thing I do for comparing assemblies is to chop one assembly into smaller pieces (e.g. 10kb), map the pieces to other other assembly to generate BAM output, then look at the BAM pileup file in a BAM viewer (e.g. Tablet).

For whole-genome align:
Maybe you could use mummer? Mummer is good to align the whole genome of species, and it also has so many other tools to modify the result. (It can show the SNP.)
For short sequence:
You could try mega X

Before you compare your genome assembly fasta, I would recommend you first do a reference guided scaffolding using the tool RagTag.
Once you are done you can compare the newly generated scaffolds (with other genomes or not) using the software BRIG. I made a video of how to use this software here
In addition you can also do a genome annotation on your scaffold using prokka. I made a video for that as well here
Good luck


text:
How to validate that BAMs have been downloaded correctly?
I currently have several hundred BAM files which were downloaded by someone else. These have remained untouched---before working with them, I would like to double-check that these BAMs have been fully downloaded.
I don't MD5 Checksums to look at.
In order to validate the BAMs, I would use something like Picard's ValidateSamFile, https://broadinstitute.github.io/picard/command-line-overview.html#ValidateSamFile
java -jar picard.jar ValidateSamFile \
      I=input.bam \
      MODE=SUMMARY

I could write a bash script, and run through these BAMs, and then read through the Picard output. That is checking whether the BAM is valid though, not whether it was downloaded fully (e.g. the BAM could be downloaded correctly, but be inherently flawed in some other way).
If the BAM is downloaded fully, isn't there some "end of file" bytes I could quickly check? Are there other more efficient methods?

samtools quickcheck is all you need. From the manual:

Quickly check that input files appear to be intact. Checks that beginning of the file contains a valid header (all formats) containing at least one target sequence and then seeks to the end of the file and checks that an end-of-file (EOF) is present and intact (BAM only).
Data in the middle of the file is not read since that would be much more time consuming, so please note that this command will not detect internal corruption, but is useful for testing that files are not truncated before performing more intensive tasks on them.
This command will exit with a non-zero exit code if any input files don't have a valid header or are missing an EOF block. Otherwise it will exit successfully (with a zero exit code).

Just a note if your cluster has older modules, or old binaries in your search path: samtools quickcheck requires samtools ≥ 1.3. For older versions you can samtools stats BAMFILE |& grep -q 'EOF marker is absent' && whatever, but this still tries to compute stats for the non-corrupted parts of the BAM so it's not the fastest solution.


text:
Is there an efficient way to check an input BAM in R?
I'm writing a function in R for an R package which takes as input a BAM.
my_func = function(my_bam){
    if (!file.exists(my_bam)){
        stop("Input a valid BAM file.")
    }

    ## do awesome stuff with BAM

}

Trying to write robust software, I would like to do more of a check whether this BAM is actually a BAM. Above, I'm just using file.exists().
I could also do something simple like "check whether the file has the extension *.bam" with
grepl('*.bam', my_bam)  ## TRUE for file.bam, FALSE for file.txt

However, is there a quick way to check that the input is actually a BAM? I would prefer to throw an error here, if the user inputs a *txt or *bedpe.
I don't think it's necessary to validate the BAM per se, whereby I believe most users rely upon Picard:
https://broadinstitute.github.io/picard/command-line-overview.html#ValidateSamFile
My first thought was Rsamtools for this operation:
https://bioconductor.org/packages/release/bioc/html/Rsamtools.html
However, it doesn't appear that there is a function here which checks the user inputs a BAM, e.g. BamFile() doesn't appear to throw an error
> library(Rsamtools)
> fake_bam = 'nonsense_string.bam'   ## a string I made up, not a real BAM
> BamFile(fake_bam)
class: BamFile
path: nonsense_string.bam
index: NA
isOpen: FALSE
yieldSize: NA
obeyQname: FALSE
asMates: FALSE
qnamePrefixEnd: NA
qnameSuffixStart: NA

How could I "check" that the user is inputting a BAM file? Perhaps pipe to samtools (which doesn't feel very stable in terms of software)?

You won’t know whether it’s a valid BAM file until you attempt reading it in full; even Devon’s method (which I recommend) only does a superficial check — i.e. it checks whether (parts of) the BAM header can be read correctly. However, this doesn’t guarantee that the file isn’t truncated or contains invalid data.
This is of course a more general problem, but it’s perhaps a bit stranger for BAMs (and other biological formats), as it’s relatively frequent to end up with technically invalid files that were produced by existing software. The only recourse is actually trying to use the file, and to handle errors appropriately.
For a very quick check™ you could also pass the file through gunzip and test whether the “magic bytes” at the beginning of the file equal 'BAM\1':
bgzf_stream = gzfile(filename, 'r')
on.exit(close(bgzf_stream))
magic = readChar(bgzf_stream, 4)
stopifnot(identical(magic, 'BAM\1'))

Take your BamFile(fake_bam) and pass it through open(). That will return a value if fake_bam is a real BAM (or CRAM or maybe SAM) file and through an error otherwise. You can then wrap all of that in a trycatch():
trycatch(open(BamFile(fake_bam)), do something, do something else)

Alternatively, ditch the trycatch() and let open() return the error if an invalid file in entered.

Just expanding upon the answers by @DevonRyan and @KonradRudolph, here's one way to put this in an R package:
Either using open(BamFile() from Devon's answer, which needs Rsamtools:
my_func = function(my_bam){

    out = tryCatch(open(Rsamtools::BamFile(my_bam)), error = function(e) e)
    if (inherits(out, "error")){
        stop("Input a valid BAM file.")
    }

    ## do awesome stuff with BAM
}

Or from Konrad's answer
my_func = function(my_bam){

    magic = readChar(gzfile(filename, 'r'), 4)
    if (!identical(magic, 'BAM\1')){
        stop("Input a valid BAM file.")
    }

    ## do awesome stuff with BAM
}

R has a function that specifically performs any(class(out) == 'error'): inherits(out, 'error'). Besides being more readable this is also potentially more efficient (though it won’t matter here).

@KonradRudolph You're right, I wasn't thinking clear. Edited.


text:
SAM format: Does the BAM "Integer or numeric array" field no longer exist? Why?
The SAM file format specification makes clear that optional fields must follow the format TAG:TYPE:VALUE
Type   Regexp matching                                       VALUE Description
A     [!-~]                                                  Printable character
i     [-+]?[0-9]+                                            Signed integer
f     [-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?                 Single-precision floating number
Z     [ !-~]*                                                Printable string, including space
H     ([0-9A-F][0-9A-F])*                                    Byte array in the Hex format10
B     [cCsSiIf](,[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)+    Integer or numeric array

However, looking at the samtools source code, I noticed this:
https://github.com/samtools/samtools/blob/develop/sam_view.c#L920-L931
switch (aux_type) {
    case 'H':
    case 'Z':
        if (kputs(bam_aux2Z(s), linebuf) < 0) return false;
        break;
    case 'i': kputw(bam_aux2i(s), linebuf); break;
    case 'I': kputuw(bam_aux2i(s), linebuf); break;
    case 'A': kputc(bam_aux2A(s), linebuf); break;
    case 'f': kputd(bam_aux2f(s), linebuf); break;
    case 'B': kputs("*** Unhandled aux type ***", linebuf); return false;
    default:  kputs("*** Unknown aux type ***", linebuf); return false;
}

Questions:
(1) I have never seen an "Integer or numeric array" of case "B". What does such an tag look like? Is there an example?
(2) Why is this currently "*** Unhandled aux type ***"? Was this the case with previous versions of samtools?

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

Thanks for the help. "which is from TMAP"---where did you find this?

Also, how command is this tag? I don't think I've seen this type of tag "in the wild" in any BAM files....though perhaps I'm very wrong.

That particular example can be found here. You don't see tags like this that frequently, I had to search for a concrete example since I don't recall seeing any locally in a few years.

Last questions: (1) I haven't seen this either. Is there any evidence that this tag format has largely been abandoned in recent years? (2) Is there a reason why the tag isn't output into FASTQ files?

That's used under the hood if you have more than 65535 CIGAR operations, but otherwise hasn't really seen wide-spread use. I guess no one ever saw the point in outputting that in fastq files. Presumably someone could make a PR to add it.

"if you have more than 65535 CIGAR operations" Could you detail this a bit more?

I'm not sure what there is to detail there. Really long noisy reads can result in a lot of CIGAR operations, which exceed the bounds of the actual CIGAR field in a BAM file. Thus the CG auxiliary tag, which is also an array at its core.

Thank you. I just needed a further clarification. I think I understand now, thank you


text:
Is there a "definitive" database for Mobile Genetic Elements?
This isn't quite my speciality, but I've been tasked with finding any evidence of Mobile Genetic Elements (like transposable elements, retrotransposons, bacteriophages, etc.) within some WGS samples.
Is there a "definitive" database used with sequences? How "species-specific" would this be?
Furthermore, are there standard algorithms used to find these MGEs? I suspect we don't have a comprehensive database of all possible MGEs, but it may be possible to query against the sequences we do indeed have (or detect any repetitive elements within the genome).

Would advise doing Denovo identification method structure-based. Then use these public databases (  RepBase, Dfam ) to look for high similarity or homology exploration, this possibility if your species, the repeats not yet discovered, or deposited in the databases.
If the species well-investigated earlier and the mobile element recognised, I assume it will be a straightforward task using species-specific search via Repeatmasker tool with the -species flag to obtain the MobileElement annotation in the genome.
Let me know if you need any data from Repbase related to your species, i have an entree.

Thanks! This is a great answer. For how long has RepBase been subscriber-only?

I think the middle of last year

Back in the Day™ RepeatMasker[1] was the go-to tool everyone used for TE annotation in eukaryotes. At the time, the RepBase database from GIRI[2] (license required) was the de facto database recommended for RepeatMasker users. I know de facto isn't necessarily definitive, but it may have been the only public database supported for some time.
It looks like more recent efforts have focused on developing the Dfam database, mirrored after the design of databases like Pfam and Rfam.

[1]Smit AFA, Hubley R, Green P (2013-2015) RepeatMasker Open-4.0, http://www.repeatmasker.org.
[2]Jurka J, Kapitonov VV, Pavlicek A, Klonowski P, Kohany O, Walichiewicz J (2005) Repbase Update, a database of eukaryotic repetitive elements. Cytogentic and Genome Research 110:462-467, doi:10.1159/000084979.

Thanks for the answer. So, if I understand correctly, the strengths of the Dfam database over RepeatMasker is that it's more extensive?


text:
How to concatenate FASTAs with respect to chromosome?
I have the following two small FASTAs, dna1.fa and dna2.fa:
# dna1.fa
>chr1
ACTGACTAGCTAGCTAACTG
>chr2
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG

# dna2.fa
>chr1
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
>chr5
CATCGATCGTACGTACGTAG

If I concatenate these two files using cat, the piped result is as follows:
$ cat dna1.fa dna2.fa > total1.fa
$ cat total1.fa
>chr1
ACTGACTAGCTAGCTAACTG
>chr2
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
>chr1
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
>chr5
CATCGATCGTACGTACGTAG

Here, as expected, these are two files which have concatenated together. It seems this is usually what is done when one concatenates several FASTAs together in bioinformatics.
Question1: Is it possible to concatenate these "based on the chromosome"? Here is the desired result:
>chr1
ACTGACTAGCTAGCTAACTG
GCATCGTAGCTAGCTACGAT
>chr2
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
ATCGATCGATCGTACGATCG
>chr5
CATCGATCGTACGTACGTAG

Note that these are exceptionally small FASTAs---usually one could be dealing with files of several GB in size.
Question2: The meta question here is, are there any standard bioinformatic analyses whereby the "concatenated by chromosome" FASTA would be more useful than the straightforward concatenated FASTA via cat alone?

Why would you do this? You're creating contiguous sequences out of (presumably) non-contiguous sequence?

Also you should not assume FASTA format is one sequence per line

What's the actual context here? Are these sequence alignments? Aligned FASTA would not follow the format of your desired output

Using the TblToFasta and FastaToTbl scripts I have posted before, combined with classic *nix utils, you can do:
$ join -a1 -a2 <(FastaToTbl dna1.fa | sort) <(FastaToTbl dna2.fa | sort) |
    sed 's/ /=/2' | TblToFasta | tr = '\n'
>chr1
ACTGACTAGCTAGCTAACTG
GCATCGTAGCTAGCTACGAT
>chr2
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
ATCGATCGATCGTACGATCG
>chr5
CATCGATCGTACGTACGTAG

Explanation
The FastaToTbl will put everything on the same line, with the sequence's ID at the beginning of the line, then a tab and then the sequence (no matter how many lines of sequence you have). So, we do this to both files, sort them and pass that as input to join, which will join its input on lines where the 1st field is the same (by default). The -a1 -a2 means "also print non-matching lines from both files" so we don't miss cases where a sequence is in only one of the two files.
Since the join will add a space between the two joined lines, we deal with that by replacing the second space on each line with a = (sed 's/ /=/2') and then pass the whole thing to TblToFasta to revert to fasta format. Finally, we use tr to replace the = with a newline so you can see where each sequence came from as you show in your example.
For longer sequences where you want to avoid sorting the input, you can do:
$ FastaToTbl dna1.fa dna2.fa |
    awk '{
            seq[$1] ? seq[$1]=seq[$1]"\n"$2 : seq[$1]=$2
         }
         END{
              for(chr in seq){
                printf ">%s\n%s\n", chr, seq[chr]
              }
         }'
>chr1
ACTGACTAGCTAGCTAACTG
GCATCGTAGCTAGCTACGAT
>chr2
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
ATCGATCGATCGTACGATCG
>chr5
CATCGATCGTACGTACGTAG

Explanation
Here, the tbl format sequences (see above) are passed to an awk script which does the work:

seq[$1] ? seq[$1]=seq[$1]"\n"$2 : seq[$1]=$2 : if the 1st field ($1, the sequence identifier) is already in the array seq, then set the value of seq[$1] to the current value (so whatever sequence has been already seen for this ID), a newline and the 2nd field ($2, the sequence). If it hasn't been seen before, then set the value of seq[$1] to the 2nd field, the sequence.
END{ for(chr in seq){ printf ">%s\n%s\n", chr, seq[chr]} : after processing everything, go through each element in the array seq and print out the data in fasta format.

However, that will mean keeping both files in memory which may be a problem for huge sequences. You can avoid that by using temp files:
$ FastaToTbl dna1.fa dna2.fa |
    awk '{
            if(!seen[$1]++){
                print ">"$1 > $1".tmp"
            }
            print $2 >> $1".tmp"
        }' && cat *.tmp > out.fa && rm *.tmp
$ cat out.fa
>chr1
ACTGACTAGCTAGCTAACTG
GCATCGTAGCTAGCTACGAT
>chr2
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
ATCGATCGATCGTACGATCG
>chr5
CATCGATCGTACGTACGTAG

Explanation
This is the same basic idea as the script above, only instead of storing the information in an array, we print it as it comes into temp files.

if(!seen[$1]++){print ">"$1 > $1".tmp"} : if this is the first time we've seen this ID, print it as a fasta header into the file $1".tmp" (e.g. chr1.tmp).
print $2 >> $1".tmp" append the current sequence to the corresponding temp file.

Once the temp files have been created, we concatenate them (cat *.tmp > out.fa) and delete them (rm *tmp).

As for why you would want to do this, an obvious example is reconstructing a reference genome from short sequences. You would want all the sequences of the same schromosome together. Of course, you would also need to assemble them in the right order, which you aren't doing here.
Alternatively, you might want to concatenate them this way so as to avoid having the extra, unneeded lines (>chrN) which would make your file larger. Or you would want to be able to retrieve all sequences with the same ID at once.
But I don't really know why you would want to do it. You're the one asking for it, after all! :)

Can you answer their question 2, the why?


text:
How to identify latent variables in single-cell RNA-Seq data
I have a single-cell RNASeq sample, in which I'd like to identify latent variables (e.g. response to stress) that I think might be affecting the clustering.
The approach I was planning to use is to reduce the dimensionality of the dataset using one algorithm among:

PCA
ICA
NMF

Then identify the most important biological processes for every component (e.g. in each Principal Component), and then finally regress out the processes that might not be relevant to the analysis or that might affect the downstream clustering.
I have a few questions:

Are there issues in this approach?
Is there a collection of gene sets (similar to MSigDB) that can be used to identify noise latent variables?
Which algorithm (PCA, ICA, NMF, or others) is more suited to this kind of analysis?

..... PCA + tSNE

@Michael why tSNE? How do you use tSNE in such case?

Its brilliant.... :-)

The larger question is how do you know you are regressing out the right stuff?

And if the so called latent variable is taking up so much of your variance... how can you be sure of the response of your interest

@StupidWolf it's a good point, however, I don't want to regress out the full Principal Component. I want to analyze which biological processes are driving the clustering and run the same analysis twice (non-regressed, and regressing only the biological processes of interest, e.g. cell cycle). My interest is more methodological (PCA, ICA and NMF work in slightly different ways, for example).

Cool.. hmm very interesting actually. I try and find the suitable scRNAseq to give it a shot. Just to clarify, yes if it is cell cycle etc we can identity it, what if it's something more ambiguous? And about regressing it out, for finding clusters etc it's not a problem. If we need to do actual inference, I would include these as covariates

I think to identify latent variables, PCA is probably not going to work. NMF might be worth trying. You might want to check out a method called consensus NMF (cNMF) (https://elifesciences.org/articles/43803) published in eLife. The author described some latent variables corresponds to cell types, while others corresponds to cell-cycle or metabolic states of cells.

Thanks! I will read it now!


text:
How to search a specific sequence in BAM files for 10X experiment
I have to search a specific sequence in a set of cells (> 1k cells) from a single-cell experiment done with 10X genomics.
As input file I have a single bam file, and 24 fastqs, therefore each file contains information about several cells. How can I search for a specific sequence and underpin which single cells have it?
I should be able to search for a specific sequence in a bam using samtools, etc. But the issue here is how to demultiplex the data from different cells to understand which cell provide the sequence.

The bam file has bam tags which say what reads belong to what cells.
https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/output/bam


text:
Is it possible to merge scRNAseq data from experiments with different design?
I have 4 different single-cell RNAseq experiments, each one representing a different sample of cell type population. I'd like to merge them to a single dataset. However, different cell types are enriched in each sample, and different samples can represent also different conditions (control and disease).
If I run a PCA on the merged dataset (without any processing of data) I find that the different sample's cells cluster separately. This could be due to biological and technical (i.e. batch effects) variation.
Is it possible to merge these different samples while retaining a good amount of information or should they be analyzed separately?

Is there some samples (controls perhaps) in both experiments? What is your goal by merging them: compare a cell line of experiment A with a cell line in experiment B?

No controls unfortunately. My goal by merging them is basically to find different sub-populations in the data and markers by exploiting a larger dataset. This data was presented to me after sequencing, and I think that it is only possible to use each experiment separately, but I'd like to be proven wrong - at least in some specific analysis.

I'm not familiar enough with the methods for scRNA, sorry. But why do you think that increasing the number of cells (of different experiments) will improve your statistical power to find new/different sub-populations? (Perhaps is specific of the methods, but I think they should already be present in each individual dataset)

I am thinking of sampling cells from different part of the same organ as a proxy for spatial resolution, for instance. The ability of merging different datasets may give the opportunity to compare the different cells populations with added spatial heterogeneity distribution. But however, this is just an example.. in my case I was wondering if merging could add more information.

How different is the experimental design? cellranger supports “reanalyse” across samples but I’m not sure this applies with different platforms generating the data.

@TomKelly they are generated with the same platform on 10X, but they are enriched differently using markers from FACS, so I think it is not possible.

@TomKelly which is the tool within cellranger to merge different samples? I was looking at https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/reanalyze but it seems to me it cannot be used to merge samples

I remember attending a seminar where scientists merged data from different types of sc experiments (i think from the same platform though). Its long ago, so my memory might not be correct on this, but I do recall that it would be a functionality planned to be implemented in Seurat - although I haven't yet had any luck finding documentation on https://satijalab.org/seurat/. In addition I think that someone at the seminar mentioned using WGCNA to map batch effects. Since my memory is not clear on this subject I don't dare post this as an answer though!

This can possibly be of interest:
https://www.ncbi.nlm.nih.gov/pubmed/29608177 - Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors

It’s right there on the cellranger manual:
#aggregate results of counts for separate samples
cellranger aggr
#analyse the combined results
cellranger reanalyze

Note: I’m not sure this applies to the above case where batch effects could be an issue. However, this is how the analysis is performed for 10X data from multiple runs in principle.
cellranger aggr
cellranger reanalyze

Adding to this answer months later: it does not work with different designs, it just runs depth normalization but it can be severely affected by batch effects.

Yes batch effects must be considered. Since this question was specifically about combining datasets, that is what was described here. This is the first step to adjust for batch effects, there are several methods to do this on an aggregate dataset and there are separate questions on this already (these assume some common cells between samples or it will overfit). I recommend not to use downsampling but aggregate raw. Depth normalisation does not remove batch effects or background noise in each sample.

I agree. As you said, methods like ComBat are not suitable for experiments with different design because they require a balanced design (some common cells between samples). Thanks.

The 3.0 release of the Seurat R package (currently in alpha) will have extensive support for combining data from different experiments.

See this paper from the Marioni group, where they propose a method for correcting batch effects between single cell sequencing experiments when each experiment contains different sub-populations:

Large-scale single-cell RNA sequencing (scRNA-seq) data sets that are
  produced in different laboratories and at different times contain
  batch effects that may compromise the integration and interpretation
  of the data. Existing scRNA-seq analysis methods incorrectly assume
  that the composition of cell populations is either known or identical
  across batches. We present a strategy for batch correction based on
  the detection of mutual nearest neighbors (MNNs) in the
  high-dimensional expression space. Our approach does not rely on
  predefined or equal population compositions across batches; instead,
  it requires only that a subset of the population be shared between
  batches. We demonstrate the superiority of our approach compared with
  existing methods by using both simulated and real scRNA-seq data sets.
  Using multiple droplet-based scRNA-seq data sets, we demonstrate that
  our MNN batch-effect-correction method can be scaled to large numbers
  of cells

Have you tried seurat3? Here is the reference:
https://satijalab.org/seurat/v3.0/pancreas_integration_label_transfer.html