text:
Why models of stochastic gene expression predict that intrinsic noise should increase as the amount of transcript decrease
I am reading Elowitz et al. (2002). When talking about intrinsic noise, i.e. noise due to microscopic events that govern which reactions occur and in what order, it is stated that:

Models of stochastic gene expression predict that intrinsic noise should increase as the amount of transcript decreases.

Noise, defined in the article as the standard deviation divided by the mean, can be divided in intrinsic (due to microscopic events during transcription) and extrinsic (due to localization, state and concentration of other cellular components, e.g. polymerases).
I think the intrinsic noise should increase at low level of transcript amount because the average (mean) amount will be lower but deviation decrease is smaller than mean decrease. I don't know, however, if this is the actual reason (and why deviation decrease should be smaller than mean decrease) or if deviation decreases with the mean and there's another more likely explanation.

Elowitz, M. B., Levine, A. J., Siggia, E. D. & Swain, P. S. Stochastic gene expression in a single cell. Science 297, 1183–1186 (2002).

Reading Kaufmann and van Oudenaarden (2007), it seems to validate the first alternative (using results from the Central Limit Theorem):

Although biochemical fluctuations influence all stages of gene expression, those involving molecules in extremely low abundance are expected from a statistical standpoint to be larger in magnitude and therefore to contribute disproportionately to the overall variation between cells.

Kaufmann, B. B. & van Oudenaarden, A. Stochastic gene expression: from single molecules to the proteome. Curr. Opin. Genet. Dev. 17, 107–112 (2007).

The inverse relationship between variance and mean is a consequence of the stochastic nature of gene expression, which can be modeled as a discrete birth-death process. One of the simplest models assumes that the gene is expressed constitutively. Then, using the chemical master equation, mRNA production can be shown to evolve as a discrete Poisson process. The steady-state distribution is Poisson. For the Poisson distribution, coefficient of variation ($\sigma^2/\mu^2$) is equal to $1/\mu$. Further, if protein is translated in bursts from mRNA (translation burst), protein noise can be written as
\begin{equation}\label{eq1}\tag{1}
\eta^2_P = \frac{\sigma_P^2}{\mu_P^2} = \frac{B}{\mu_P},
\end{equation}
where $\eta_P^2, \sigma_P^2,\mu_P,$ and B are the coefficient of variation, variance, mean and burst size for protein. From \eqref{eq1}, the inverse relationship is evident.
Gene expression models can be made more detailed by assuming different states for the promoter. The promoter switches between these states, and each of these states can have different rates of transcription. For these models, $\eta_P^2$ can again be shown to follow an inverse relationship with $\mu_P^2$.
Reference:
Paulsson, Johan. "Models of stochastic gene expression." Physics of life reviews 2.2 (2005): 157-175.


text:
Demultiplexing and preprocessing for custom single-nucleus Drop-seq data
I am trying to reproduce the preprocessing of paired-end sequencing reads in Lake et al. (ref. 1).

First, paired-end reads were filtered out if read 1 had more than four non-T bases in the last ten bases (to remove all non-poly(T)-captured contaminated reads), or had one or more bases with a poor quality score (<10). Cell barcode and UMI information were then inferred from the first 12 bases and the next 8 bases of read 1, respectively. The right mate of each read pair was trimmed to remove any portion of the SMART adaptor sequence or large stretches of poly(A) tails (6 consecutive bp or larger).

I was able to filter out reads as described in the first sentence.
However, I am unsure on which tool or procedure to use to demultiplex based on cell barcode and UMI, and to trim the right mate. Do you know what I can use or something I can read to understand how to do it?

Lake, B. B. et al. Integrative single-cell analysis of transcriptional and epigenetic states in the human adult brain. Nature Biotechnology 36, 70–80 (2017).

I didn't notice the sentence:

Paired-end sequencing reads were processed largely as described (http://mccarrolllab.com/wp-content/uploads/2016/03/Drop-seqAlignmentCookbookv1.2Jan2016.pdf ), with additional correction steps.

This should provide instructions on how to lead the preprocessing and mapping of data.

source: Stack Exchange

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

sort: votes


text:
Which tools for differential expression analysis in scRNA-Seq?
I am starting to run analysis for differential expression in scRNA-Seq. Which tools are available for this kind of analysis?
Can tools for bulk RNA-Seq like DESeq be used for scRNA-Seq?

I recommend this page to learn about other tools

@Llopis thanks. I'd recommend also this curated repository for Single-Cell analysis: https://github.com/seandavi/awesome-single-cell

Yes, and there's also the online course, where it says that "Several computational analysis methods from bulk RNA-seq can be used"

According to Jaakkola et al. (ref. 1):

Based on our comparisons, DESeq and Limma without any modifications are not suitable for scRNA-seq data analysis, and yet, they have performed well in the context of bulk RNA-seq data.

They also provided a summary of the comparison between 5 tools (2 designed for Single-Cell, 2 designed for bulk, 1 general):

In my case I ended up using MAST (ref. 2, github).

Jaakkola, M. K., Seyednasrollah, F., Mehmood, A. & Elo, L. L. Comparison of methods to detect differentially expressed genes between single-cell populations. Briefings in Bioinformatics 18, 735–743 (2017).
Finak, G. et al. MAST: A flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data. Genome Biol. 16, 1183 (2015).

In a paper of Soneson & Robinson (2018) in Nature Methods (doi:10.1038/nmeth.4612), the authors concluded that "bulk" methods do not perform worse than methods developed specifically for scRNA-seq.

@b.nota Thank you! Can you edit my question to reflect also your comment?

license: 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:
How to apply function on bedfile based on a reference bedfile
I have two bed-like files:
File1
chr1     1     10    .    1    Sample1
chr1    20    180    .    1    Sample1
chr1     1      5    .    1    Sample2
chr1    70    100    .    1    Sample3

File2 (reference)
chr1     1     5     .    0
chr1     5    20     .    1
chr1    50   120     .   -1

I'd like to find a way to apply a function (in my case a subtraction) between the first file and the second, using the second as a reference, and producing a results file containing also the indication of the sample name.
For example, in my case the results file would be:
chr1     1      5    .    1    Sample1
chr1     5     10    .    0    Sample1
chr1    20     50    .    1    Sample1
chr1    50    120    .    2    Sample1
chr1   120    180    .    1    Sample1
chr1     1      5    .    1    Sample2
chr1    70    100    .    2    Sample3

I was thinking about implementing it in pandas, but it would be long to do and probably not so efficient. I looked for bedmap but it seems to me that this operation can be carried on (with the sum operation) only within a single file during aggregation, not between two different bed files. So maybe it is not the right tool to use in this case.
Also, it would be important in the final file to retain the sample name.

You could munge your first file into a BED5 file:
$ awk -vOFS="\t" '{ print $1, $2, $3, $6, $5 }' /tmp/f1.bed | sort-bed - > /tmp/f1.bed5

Sort the second file, if necessary:
$ sort-bed /tmp/f2.unknownSort.bed5 > /tmp/f2.bed5

Take the union of the two sorted files:
$ bedops -u /tmp/f1.bed5 /tmp/f2.bed5 > /tmp/union.bed

Here's what the union looks like:
chr1    1       5       .       0
chr1    1       5       Sample2 1
chr1    1       10      Sample1 1
chr1    5       20      .       1
chr1    20      180     Sample1 1
chr1    50      120     .       -1
chr1    70      100     Sample3 1

We have one name scheme for reference element IDs (.) and a second name scheme for sample IDs (Sample*).
We will exploit this difference in name schemes when we apply operations in a certain order. For this example, we will subtract a reference element score from a sample element score.
We can partition the union, mapping the disjoint elements back to the unioned set:
$ bedops --partition /tmp/union.bed | bedmap --echo --echo-map-id --echo-map-score --delim '\t' - /tmp/union.bed
chr1    1   5   .;Sample2;Sample1   0.000000;1.000000;1.000000
chr1    5   10  Sample1;.   1.000000;1.000000
chr1    10  20  .   1.000000
chr1    20  50  Sample1 1.000000
chr1    50  70  Sample1;.   1.000000;-1.000000
chr1    70  100 Sample1;.;Sample3   1.000000;-1.000000;1.000000
chr1    100 120 Sample1;.   1.000000;-1.000000
chr1    120 180 Sample1 1.000000

Each of the fourth and fifth columns are semi-colon-delimited strings. Each element is ordered by the order in which it is discovered in a bedmap operation.
Both columns have the same ordering or same one-to-one correspondence of ID to signal score. This is very useful, because we can turn these into a pair of ordered arrays and apply some logic, depending on the answer we want back.
For example, here is some awk logic to parse those two strings and subtract a reference element score from a sample element score:
$ bedops --partition /tmp/union.bed \
    | bedmap --echo --echo-map-id --echo-map-score --delim '\t' - /tmp/union.bed \
    | awk -vFS="\t" -vOFS="\t" \
        '{ \
             i=split($4,a,";"); \
             j=split($5,b,";"); \
             if ((i==1) && (a[i]==".")) { pass; } \
             else { \
                 for (k=1; k<=i; k++) { \
                     if (a[k]==".") { refIdx=k; } \
                 } \
                 for (k=1; k<=i; k++) { \
                     if (k!=refIdx) { print $0, ".", b[k]-b[refIdx], a[k]; } \
                 } \
             } \
         }' \
    | cut -f1-3,6-8 \
    | sort-bed - \
    > answer.bed

The file answer.bed will look like this:
chr1    1   5   .   1   Sample1
chr1    1   5   .   1   Sample2
chr1    5   10  .   0   Sample1
chr1    20  50  .   1   Sample1
chr1    50  70  .   2   Sample1
chr1    70  100 .   2   Sample1
chr1    70  100 .   2   Sample3
chr1    100 120 .   2   Sample1
chr1    120 180 .   1   Sample1

This partitioning of intervals is different from the one in your question. Perhaps this is due to how you partitioned your set. In any case, perhaps this would give you an idea of how to apply custom operations over a disjoint partitioning.
The main idea is that this part of the awk script finds the score associated with the reference file ("File2"):
for (k=1; k<=i; k++) {
    if (a[k]==".") { refIdx=k; }
}

Any index other than refIdx corresponds to a Sample* element. That means we can do operations between Sample* scores and the specific reference score.
We also leave out mappings where a reference score does not associate with any Sample* element:
if ((i==1) && (a[i]==".")) { pass; }

Thus, all the operations are always between a pairing of a sample and the mapped reference element.
The cut and sort-bed statements at the end just clean up the result so that you get something that looks like what you are expecting at the end, which is also sorted correctly and ready to use for downstream set operations.


text:
TagReadWithGene missing when using latest version of Drop-seq_tools
I am using Drop-seq_tools to analyze scRNA-seq in a similar way of this paper. I am using the same data (GSE97930) to get used to Drop-seq_tools.
On the GEO repo the authors provide the relevant code (GSE97930_README.txt). However, one tool used in the script is missing: TagReadWithGene.
The latest version of Drop-seq_tools does not have such method (as it should be provided by this library), and there is no version history for the whole package, so I cannot download previous versions.
The best matches for the missing method are TagReadWithGeneExon and TagReadWithInterval.
In the methods of the cited paper they say:

Reads that mapped to intronic or exonic regions of genes as per the GENCODE gene annotation were recorded.

I think that in this case the tool to use should be TagReadWithInterval since the other one should consider only exons.
I cannot find any other documentation about the tools. Do you know/already encountered this problem?

As per this response from Drop-seq_tools team, at the moment intronic regions cannot be annotated and TagReadWithGeneExon should be used instead.
However, in the version of the public repository:

TagReadWithGeneExonFunction is the old code, for people who want the more simple tagging methodology and rely on it for 3rd party tools, and TagReadWithGeneFunction, which for your purposes replaces TagReadWithGeneExon. [...]  You’ll see additional options to extract expression for CODING, UTR, INTRONIC reads [...]


text:
Does read sequence identifier need to be unique for proper alignment?
I am preprocessing a set of scTHSSeq reads in which each read sequence identifier has been replaced with the cell barcode, e.g.:
@CGAGGCTG_TATCCTCT_AGTCCA
...
...
...
@CGAGGCTG_TATCCTCT_GGTACG
...

The read sequence identifier has been replaced with the cell barcode to be able to demultiplex alignments after bwa into a different bamfile for each cell.
I ran the alignment using bwa mem and I received no errors. However I would like to know if there are caveats or reasons why this approach may lead to misleading results.
I will then proceed to remove PCR duplicates with samtools rmdup executed in each cell's bamfile and then I'd run peak calling using spp.

What tool did you use to replace the read names? Can you modify that to split up the reads into separate files or encode the cell barcode in the read description instead of the read name?
Typically is it assumed by many tools that the FASTQ read names are unique. you could have issues with BWA and split alignments but I don't know if that pertains to your data.  I think the best you can do is the sanity check the results at each step of the process.
In the fastq format (https://en.wikipedia.org/wiki/FASTQ_format) on the first line denoted by the '@' there can be two values separated by a space.  The first is the name and the second is the description.
So an example would be:
@Original_read_name CGAGGCTG_TATCCTCT_AGTCCA
ATCGATCATGATGACA
+
@&%*@(%@%*($!##!

In this iteration I did it using awk, do you mean incorporating it in the first line after the sequence identifier? I can certainly modify the script.

I mean the one defined as entry.comment here?

I edited my answer to address your question

I noticed now that BWA does not preserve read description, I should probably append cell barcode to read name or this may introduce problems?

The classic way to do this would be to append the barcodes to the read name such that the read name looks like: @Original_read_name_CGAGGCTG_TATCCTCT_AGTCCA, If the barcodes are extacted from the read sequences themselves, this processing can be achevied using umi_tools extract


text:
How is PCR duplication rate computed in scATAC-seq?
Reading Cusanovich et al. (2015) I encountered the sentence:

We mixed pairs of cell lines (HEK293T or HL-60 versus GM12878), performed combinatorial cellular indexing, and sequenced the resulting libraries to saturation (65% duplicate rate).

How is PCR duplication rate computed?

Cusanovich, D. A. et al. Multiplex single-cell profiling of chromatin accessibility by combinatorial cellular indexing. Science 348, 910–914 (2015).

The answer is in the supplementary methods page 3.

Finally, all fragments in the same library with duplicate start and
  end coordinates were removed using Picard


text:
How to detect mismatch before mapping in RNA-Seq data
In the methods of this paper, the authors say:

Reads were then filtered based on quality score in the UMI region. Any read with >1 low-quality base (phred <=10) were discarded. Reads with more than one mismatch in any of the three 8 nt cell barcodes were also discarded. The cDNA reads (Read 1) were then mapped to either a combined mm10- hg19 genome or the mm10 genome using STAR (44).

I understand how to filter based on mismatch after mapping but in this case it is done before mapping. Has this to be interpreted as "Reads with low quality in more than one base in any of the three 8 nt cell barcodes were also discarded"?

At each step in the SLiTing and pooling the barcodes added are not random, but are taken  from a predefined list of possible barcodes.
When the reads arrive, first you filter by reads that have >1 low quality bases in the UMI region (which is truely random).
Then you compare the sequences in the 3 cell-barcode regions against the list of allowed barcodes for each of those positions. If a read is not within an edit distance of one from one of those barcodes, it is discarded.
Finally the reads that pass these filters can be aligned.

It should not be interpreted as "Reads with low quality in more than one base in any of the three 8 nt cell barcodes were also discarded". They only care if the cell barcodes match what are valid cell barcodes within 1 mismatch, the individual qualities are ignored. This sort of filtering turns out to be pretty common in practice, since you have to write a maximum edit distance into whatever you use to demultiplex by cells (or labeling each read, as was likely done in this case).

Cell barcodes =/= sequence reads.  The default settings on Illumina's bcl2fastq work exactly the same way; they allow a read to count for a sample if there is a single mismatch in the index, but not more.


text:
Definition of "proband" in CNV analysis
I am reading this paper regarding de novo CNV analysis in Schizophrenia:

Accordingly, such CNVs must frequently occur as recurrent de novo mutations. In a sample of 662 schizophrenia proband–parent trios, we found that rare de novo CNV mutations were significantly more frequent in cases (5.1% all cases, 5.5% family history negative) compared with 2.2% among 2623 controls, confirming the involvement of de novo CNVs in the pathogenesis of schizophrenia.

I assumed that proband stands for the sample being tested (child of the parents), however, I couldn't find an explicit definition for it. Is it just the sample being tested, or does it imply additional details?

answer from @maximilian-press, converted from comments:

https://www.genome.gov/genetics-glossary/Proband
https://www.merriam-webster.com/dictionary/proband

A proband is an individual affected with a disorder who is the first subject in a study (as of a genetic character in a family lineage) who brings a genetic disorder to the attention of the medical community.


text:
Definition of Copy Number Abundance
Reading from this paper:

The remaining 9 new copy-number alterations showed altered coverage patterns on whole-genome sequencing but could not be confirmed by orthogonal methods because of their small size, low abundance, or both (Fig. 2C, Fig. S2C, and Table S5)

In the supplementary file, CNA abundance as percentage is defined as:
$$Abundance = \frac{2^{L2R} - 1}{\frac{CN}{2} - 1}*100$$
where L2R is the log$_{2}$ normalized coverage ratio vs. a panel of normals and CN is the estimated copy number for the event.
I haven't found a better definition/explanation behind "CNA abundance", is it defined elsewhere, or is it a common term in CNA analysis?

My guess by analyzing the definition is that CNA abundance is a proxy for the percentage of (tumor) cells carrying the alteration.


text:
On the same strand, for the same gene, can exons be overlapping?
I want to get a set of exon regions for each protein coding gene.
I extracted a set of relevant information (chromosome, start, end, gene ID, gene name, gene type, exon number and exon ID) from a GTF using:
awk '$3 == "exon"' annotation.gtf \
  | gtf2bed \
  | cut -f1-4,10 \
  | awk 'match($0, /gene_name \"[a-zA-Z0-9\_\-.\/]*\"/, gene_name)
    && match($0, /gene_type \"[^"]*\"/, gene_type)
    && match($0, /exon_number [0-9]*/, exon_number)
    && match($0, /exon_id \"[^"]*\"/, exon_id) {
      print $1,$2,$3,$4,gene_name[0],gene_type[0],exon_number[0],exon_id[0];
    }' \
  | tr -d '"' \
  | tr ' ' $'\t' \
  | cut -f1-4,6,8,10,12 \
  | grep 'protein_coding' \
  > annotation.exons.bed

However, during a sanity check, I noticed that for some genes there are some exons with overlapping coordinates, same gene information and exon number, but different exon ID, e.g.:
chr2    9423650 9424052 ENSG00000119203.14      CPSF3   protein_coding  1       ENSE00001829522.1
chr2    9423653 9423823 ENSG00000119203.14      CPSF3   protein_coding  1       ENSE00001830696.2
chr2    9423883 9424960 ENSG00000119203.14      CPSF3   protein_coding  1       ENSE00001932298.1

I did not directly work with exons in the past, but I am aware of alternative splicing.
My understanding was that exons can be skipped during alternative splicing, but I assumed that exon number was unique for bases, within the same gene (on the same strand).
In the example above, the first 2 listed exons have the same exon number, overlapping (but not exactly the same) bases, and different exon ID. Is this due to some additional information that I can retrieve from the GTF, i.e., from:
chr2    9423650 9424052 ENSG00000119203.14      .       +       HAVANA  exon    .       gene_id "ENSG00000119203.14"; transcript_id "ENST00000475482.5"; gene_type "protein_coding"; gene_name "CPSF3"; transcript_type "protein_coding"; transcript_name "CPSF3-203"; exon_number 1; exon_id "ENSE00001829522.1"; level 2; protein_id "ENSP00000419744.1"; transcript_support_level "4"; hgnc_id "HGNC:2326"; tag "mRNA_end_NF"; tag "cds_end_NF"; havana_gene "OTTHUMG00000090415.5"; havana_transcript "OTTHUMT00000353331.3";
chr2    9423653 9423823 ENSG00000119203.14      .       +       HAVANA  exon    .       gene_id "ENSG00000119203.14"; transcript_id "ENST00000238112.8"; gene_type "protein_coding"; gene_name "CPSF3"; transcript_type "protein_coding"; transcript_name "CPSF3-201"; exon_number 1; exon_id "ENSE00001830696.2"; level 2; protein_id "ENSP00000238112.3"; transcript_support_level "1"; hgnc_id "HGNC:2326"; tag "basic"; tag "MANE_Select"; tag "appris_principal_1"; tag "CCDS"; ccdsid "CCDS1664.1"; havana_gene "OTTHUMG00000090415.5"; havana_transcript "OTTHUMT00000206843.2";

?
Is this because the different exons belong to different transcripts?
However, shouldn't different exons still be non-overlapping one against each other?

I'm not quite sure what you're asking, but pretty much anything can happen in terms of exons and transcripts:

Two different genes encoded on different strands can share the same genomic region as an exon

Two different genes encoded on the same strand can share the same genomic region as an exon (including identical exon boundaries)

A gene can have multiple transcript isoforms that share the same genomic region on the same frame

A gene can have multiple isoforms that share the same genomic region in different frames

Genes can have variable regions which are modified during somatic development (e.g. T-cell and B-cell receptors)

Gene transcripts can be spliced together in different orders in different isoforms

Exons from one transcript isoform can overlap with exons from another isoform (but not necessarily cover identical regions)

Stop codons can occasionally be "read-through", resulting in a longer protein/peptide than the expected translation

As an example of different isoforms within the same gene exhibiting different translation mechanisms, here are graphical representations of annotated transcript isoforms for Gapdh on the mouse GRCm38 assembly, created in JBrowse. I wouldn't put complete trust in any of these (especially not the ones that have no defined intronic regions), but it's a useful demonstration of the types of things that can happen. Compare, Gapdh-207, where an untranslated region starts halfway through an exon, and Gapdh-209, where the translated protein starts halfway through an untranslated exon in Gapdh-207:

If there is a situation where it seems like exons from the same isoform have overlapping regions, it might be worth looking elsewhere in the genome to see if there's a transcribed gene copy somewhere. That, at least, doesn't seem to be something I'd expect to happen through the normal translation process.
Also, don't expect annotations to be complete. Even on something as [presumably] well-studied as the mouse mitochondrial genome, we can identify clear transcript start/stop points from long polyA cDNA sequences that are completely absent in the genome annotations. This is most obvious for the completely unannotated reverse strand (e.g. transcript starting at around 16,200bp), but also happens on the forward strand as well (e.g. coverage steps within mt-Nd5):

Thank you, I believe my confusion was due to the fact that I was not aware of the exon behavior from your two last points. This I think might explain what it is happening in my case, where the exon 1 from the first transcript shares some bases (but not exactly the same interval) with exon 1 from the second transcript. I used to think that transcripts (on the same strand) might have different exons, but the set of exons from which they could select their subset was composed by non-overlapping exons.

May I know which software or tool you are using to show Gadph exons in the images your provided here

I've updated my answer. I used JBrowse (v1 to be specific, but v2 is the most current version).


text:
Potential side effects of replacing read group tags in BAM file
I have a set of BAM files where the read group tags have some (default?) values, i.e.:
@RG ID:RG0  LB:LB0  PU:PU0  SM:SM0

This creates issues in my downstream analyses, where multiple BAM files with the same SM tag are used.
Samtools provides a command to replace the read group tag. However, I am not sure if there are possible side-effects that I should be aware of, and I might need to remap the BAM file after this replacement. I want to replace ID and SM with the sample id, which is unique for each BAM file.
Do I need to remap and/or run additional steps, or replacing the RG tags should be sufficient to update the BAM files in a consistent way?

potential (not exactly?) duplicate of https://bioinformatics.stackexchange.com/q/908/57. Let me know if those answers fix your problem too.

Thanks @KamilSJaron, as indicated above, my question was mainly concerned on potential side effects of using samtools to modify the header - and if remapping was needed. Therefore, even if the question you linked might be interesting, it is not exactly a duplicate of that. Thanks for your help, though.

You do not need to remap the files, replacing the read group information with samtools is sufficient to deal with this. When you update your pipeline, have it use the sample name to construct the read group information. That way you won't run into this problem again.

I assume you are using GATK software tools for downstream analysis. If that's the case, then you will most likely encounter issues due to missing information on read groups. I believe that GATK can't differentiate between read groups if they are all similar.
In all cases, you can find more information about SAM/BAM and related specifications on the 1.3 The header section as they describe what is mandatory.
So, if you have access to FastQ files, you can add them directly during the alignment process with bwa mem:

bwa mem

-R STR    Complete read group header line. ’\t’ can be used in STR and will be converted to a TAB in the output SAM. The read group ID will
be attached to every read in the output. An example is
’@RG\tID:foo\tSM:bar’. [null]

An example command line would be:
$bwa mem \
-M \                 // Mark shorter split hits as secondary (for Picard compatibility).
-t "$cpus" \         // number of threads
"$genome" \          // genome reference
"$fastq_r1" \        // R1
"$fastq_r2" \        // provide this line if you have R2
-R "@RG\tID:1\tLB:lib1\tPL:illumina\tSM:$sampleName\tPU:runBarcode" \
> "$sampleName".sam  // Output

bwa2 mem

Bwa-mem2 is the next version of the bwa-mem algorithm in bwa. It
produces alignment identical to bwa and is ~1.3-3.1x faster depending
on the use-case, dataset and the running machine.

The command line is exactly similar to the command line above, except that you would be using bwa2.

Samtools

As you have mentioned, you can completely overwrite (default mode) the read groups already existing by doing:
samtools addreplacerg \
-r "@RG\tID:1\tLB:lib1\tPL:illumina\tSM:$sampleName\tPU:runBarcode" \
-m overwrite_all \         // mode to replace all
-@ 5 \                     // number of cpus
-O b \                     // output bam
-o "$sampleName"_RG.bam \
"$sampleName".bam

Picard AddOrReplaceReadGroups

java -jar $picard AddOrReplaceReadGroups \     // Also available in GATK tools
I="$sampleName".bam \
OUTPUT="$sampleName"_RG.bam \
RGID=1 \
RGLB=library \
RGPL=illumina \
RGPU=runBarcode \
RGSM=${sampleName} \
SORT_ORDER=coordinate

Thank you for your detailed answer. My question was actually more limited to: do I need to remap or update the bam files after having modified the read group with samtools? Therefore, I decided to use the other answer. I will however award the bounty to both once it is possible for me to do it.

@gc5 you can only award a bounty to a single answer.

Hi @terdon, yes. Although I meant to award a second bounty for the help. I will be awarding it in 23 hours when it is possible for me to do it.


text:
Definition of genotype in demuxlet
I am reading the Online methods of the demuxlet paper.
The genotype $g$ is taken from the set $\{0, 1, 2\}$, defined as "the true genotype of the sample corresponding to $c$-th droplet at $v$-th variant".
What is the mapping between these values and the corresponding genotype?
Each variant is modeled using a binomial with $n=2,p=\frac{g}{2}$, so I would assume that 1 is heterozygous, 0 is homozygous REF and 2 homozygous ALT. Is this assumption correct?

The genotypes $\{0,1,2\}$ are mappings to the genotypes AA (homref), AB (heterozygous), BB (homalt) (from supplementary resources). This is also mapped to the number of alterations (e.g., 0 -> 0 alterations -> homref).


text:
Meaning of order of alleles in VCF GT field for unphased genotypes
Order of alleles in genotype field (on VCFs) can be used to infer haplotypes
However, what about unphased genotypes? I have a VCF file where I only have unphased genotypes:
0/1
1/0
1/1
1/2
2/1

However, I have, for example, both 0/1 and 1/0. Does the order of alleles have some different meaning with unphased genotypes, or genotypes with same alleles (e.g., 0/1 and 1/0) can be considered synonymous?

The fact that they are unphased necessarily means there is no order to them. However, the convention is usually to place the smaller value first, like 0/1 or 1/2 etc. But in reality it doesn't matter.

In the case of unphased VCF 0/1 and 1/0 are exactly equivalent and the order has no meaning.
Throughout the VCF specification, it is always encouraged to use numerical precedence in order to order the genotypes, e.g. 0/0, 0/1, 1/1 would be the preferred representations for diploid bi-allelic genotypes.
However it is not strictly enforced in the VCF specification as for now and you could have 1/0, 2/0, etc, and bcftools will not complain parsing them. The presence of two representation in the same VCF might indicate that the data comes from multiple variant callers or has been through manual editing, since for example variant callers will always use one representation only (numerical precedence usually)
Also, the order does impact how you read some of the sample information fields, for example the AD (allelic depth), for which the values are read following the order present in the GT field.
I had a look at the VCF 4.2 and 4.3 spec, but indeed is not so precise there on these details.


text:
Alternatives to imputation for single-cell RNA-Seq analysis to increase statistical power of DE analysis
I have a small scRNA-Seq dataset (n = 357, inhibitory neurons). This set of cells is split almost evenly between two conditions (Case and Control). I would like to test for differential expression with MAST [1]. However, with such a small number of cells and high heterogeneity of different inhibitory neuron cell types, I am able to get only a single gene as differentially expressed (False Discovery Rate 0.05, Log Fold Change 0.15).
The best solution would be to enrich this subset of neurons in a new sample, however, are there other approaches that can be used to extract more information from this dataset? I'm not considering imputation at the moment since it seems to produce a large number of false positives [2].

Finak, Greg, et al. "MAST: a flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data." Genome biology 16.1 (2015): 1-13.
Andrews, Tallulah S., and Martin Hemberg. "False signals induced by single-cell imputation." F1000Research 7 (2018).

Maybe your Case and Control cells are simply not that different transcriptionally? Have you examined their separation with t-SNE or similar?

@Chris_Rands good point, I think your suggestion explains the small log-Fold-Change. However, the Case is composed of cells with a specific heterozygous deletion, therefore I would expect to recover as differential at least the deleted genes. This is indeed what is validated in bulk RNASeq. Differently, in single-cell data, we noticed a drop in differential expression almost perfectly correlated with the decreasing number of cells. This is expected, but before running an additional sample with enriched cell types we wanted to look for existing techniques to recover part of the signal.

@Chris_Rands to reply directly to your question, they are not separated on visualization projections such as UMAP or t-SNE. We don't expect, however, a very big separation since we expect the processes involved in the deletion to be more subtle than processes associated with cell-type identity. To conclude, when we have bigger clusters of cells, we are able to recover (positive control) the genes deleted in the Case.


text:
Fast filtering of intervals not falling within a certain distance from known genes
I would like to filter a bed file with intervals, ie. in the format of:
chr1    13800   14301
chr1    15500   16001
chr1    19400   19901
chr1    22800   23301

In particular, I want to filter out intervals that fall far from known genes, let's say > 1KB.
I was thinking of querying each interval (in parallel) and check through a known set of coordinates from a gtf, e.g. from Homo_sapiens.GRCh38.93.gtf by filtering for gene coordinates:
1       havana  gene    168317497       168376876       .       +       .       gene_id "ENSG00000227722"; gene_version "2"; gene_name "AL022100.1"; gene_source "havana"; gene_biotype "transcribed_processed_pseudogene";
1       mirbase gene    168375524       168375621       .       +       .       gene_id "ENSG00000207974"; gene_version "1"; gene_name "MIR557"; gene_source "mirbase"; gene_biotype "miRNA";
1       havana  gene    168400829       168495685       .       -       .       gene_id "ENSG00000228697"; gene_version "2"; gene_name "AL023755.1"; gene_source "havana"; gene_biotype "lincRNA";
1       havana  gene    168401483       168407274       .       -       .       gene_id "ENSG00000283255"; gene_version "1"; gene_name "AL022100.2"; gene_source "havana"; gene_biotype "lincRNA";
1       havana  gene    168449672       168451211       .       -       .       gene_id "ENSG00000215833"; gene_version "3"; gene_name "QRSL1P1"; gene_source "havana"; gene_biotype "processed_pseudogene";

However, since the bedfile is quite big I am looking for the fastest solution since checking each interval may be quite expensive.

In one line, using bedtools
zcat  Homo_sapiens.GRCh38.93.gtf.gz \
 | awk '$3=="gene"' \
 | bedtools slop -b 10000 -g contigs.tsv -i - \
 | bedtools intersect -u -a intervals.bed -b -

This first takes the genes and filters them for the third column being gene. This is important because all GTF lines contain the word gene as gene_id is a manditory attribute for GTF.
Next it passes it to slop, which adds 1 kb either side (ensuring not to go below 0 or off the end of the chr: this isn't that important because we are never outputting these regions anyway, but bedtools requires it). We tell bedtools to work on the standard in by using -i -.
Finally we intersect these gene regions with the intervals.bed. -u says output the line from -a if it overlaps with at least one entry in -b and passing -b - means read file b from the standard in.
This will only work if your GTF contains gene lines, which the GENCODE/Ensembl ones do. However, it is interesting to note that these GTFs are not actaully valid GTF files, as to be valid a GTF file must have both a gene_id and transcript_id attribute, and the gene lines in the GENCODE GTFs to not have this.
If your GTF does not contain gene lines, but can create them by merging all the exons in all the transcripts and filling in the gaps using gtf2gtf from cgat.
 cgat gtf2gtf --method=merge-transcripts --stdin Homo_sapiens.GRCh38.93.gtf.gz \
 | bedtools slop -b 10000 -g contigs.tsv -i - \
 | bedtools intersect -u -a intervals.bed -b -

Here's a way to use BEDOPS, which was designed to work fast by using sorted input. Other tools now use sorting to accomplish similar performance benefits.

Convert GTF annotations to a sorted BED file of genes:
$ awk '($3=="gene")' annotations.gtf | gtf2bed - > genes.bed

Sort your intervals, if unsorted:
$ sort-bed intervals.unsorted.bed > intervals.bed

Filter intervals against 1kb-padded genes:
$ bedops --merge genes.bed | bedops --range 1000 --everything - | bedops --not-element-of 1 intervals.bed - > intervals.filtered.bed

We first --merge genes to simplify dealing with overlapping genes and further speed downstream filtering.
The file intervals.filtered.bed will contain intervals from intervals.bed which lie outside of a 1kb-window around merged gene space.

To run all of this as a one-liner (assuming use of the bash shell):
$ bedops --merge <(awk '($3=="gene")' annotations.gtf | gtf2bed -) | bedops --range 1000 --everything - | bedops --not-element-of 1 <(sort-bed intervals.unsorted.bed) - > intervals.filtered.bed

This creates the same output result, but avoids creation of intermediate files and so should run even faster.

Note: BEDOPS is agnostic about chromosome sizes — this toolkit works with any genome assembly. However, using --range takes care of ranges that would have a negative start position, setting their start coordinates to 0. Bounds aren't really relevant here, in any case, because the bedops --not-element-of step queries input intervals, which should already be within chromosome bounds, unless they were generated incorrectly.

I think there's a missing quote on the first awk call

When executing 1 I get Error: Potentially missing gene or transcript ID from GTF attributes (malformed GTF at line [1]?)

https://www.biostars.org/p/206342/#206350

First I prepared a bed file in which the gene intervals are augmented by 1KB before and after the gene start and end coordinates. Then I intersected this bed file with my original one with the option -wa, therefore retaining only the intervals in my original bed file that intersect with the bedfile produced from the original gtf, filtered by gene regions. The -u flag is also used to print only one line per overlap between the two files, in case an interval in the first overlaps with multiple intervals in the second.
# download information about chromosomes
wget http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.chrom.sizes ./

# Original GTF has a different notation than UCSC
cat Homo_sapiens.GRCh38.93.gtf | \
  sed -e $'s/^./chr&/' -e $'s/^chrMT/chrM/' | \
  awk '$3=="gene"' | \
  bedtools slop -b 1000 -g hg38.chrom.sizes -i - | \
  sed -e $'s/\tgene_id.*gene_name/\t/' -e $'s/gene_source.*//' | \
  tr -d '"; ' | \
  awk -v OFS=$'\t' '{print $1, $4, $5, $9}' > genes_1kb.bed

bedtools intersect -wa -u -a intervals.bed -b genes_1kb.bed > intervals.1kb.bed

note that this will filter out intervals that overlap a gene, or 1kp upstream, but will leave features that are within 1kbp downstream. And it will only work if the GTF has gene lines, which many GTF files don't have.

@IanSudbery, why this will filter out intervals that overlap a gne or 1kb upstream?

Oops I noticed now I forgot to add 1kb to the end. Thanks! Does it solve the problem?

@gc5 your awk command could have been replaced with bedtools slop ref which has the added advantage of making sure the intervals do not extend more than the chromosome length. It also reduces the chances of zero or one base coordinate errors (tried and test tool).

@fridaymeetssunday amazing! I didn't know it, thanks


text:
Strategies to disentangle the effect of two genes on phenotype
I am analyzing the effect of one putative tumor suppressor (and its loss) on a specific cancer phenotype. This putative tumor suppressor is known in literature to have an effect on the growth of tumors in mouse models.
However, in the same region of this gene there's a second well-known tumor suppressor - actually a very powerful one. By looking at the data from cBio portal, there's a very strong co-occurrence of loss on both tumor suppressors in a pan-cancer analysis, that is they are almost always lost together.
Therefore, I would like to disentangle the effect of the two genes, and to analyze only the effect of my putative candidate on the phenotype. The simplest approach is to use only the samples that have a single loss (only the loss of my putative candidate, but not the loss of the confounding tumor suppressor). However, only ~10% of all samples shows such behavior, and this 10% is distributed across all cancer types, decreasing the statistical power of the analysis - since I would like to test each cancer type separately to avoid mixing incompatible data.
Which kind of strategy/analysis can be used in this case?

How do you model the effect of the genes? I would use all the data (or cancer specific data) and indicate which sample have lost the gene of interest, the other or both genes. That way you could disentangle the effect of each gene and the interaction of both

Maybe you could fit a mixed model to predict phenotype, using the two genes as fixed effects and cancer type as a random effect, and ask the model for the separate effects of each gene, and the interaction between them.

Or even more simply using generalized linear modeling (GLM), where the regression weights are the genes and the accuracy of the model is measured in the residual. You also Q-Q plots to assess the fit of the model to the data.
If the residual is high and Q-Q "wobbly" transformations works wonders.

Machine learning such as a neural network is a very cool way to do this calculation, but IMO you need Python - although people on this board do this in R. The concept is very easy you are catagorising (keep in mind I'm not an oncologist) gene vs. cancer resistance and you want neural network, linear-based regression etc.. to make the call (you will see it visually). You then use the weights to assess the relative contribution of each gene.
The 10% issue is formally assessed, its not hard to obtain the power because its a common problem. The formal way to assess this, from memory, is you assign a "null model" and the accuracy of a ML model needs to a lot higher than the null model. I forget the details, I've never been in these situations. There's a specific technical name (which I should look up).

The explicit GLM / GLMM (non-ML) is the easiest approach however. If you use a Mac there is a great (relatively cheap) package for GLM, but any good software package will do. I assume GLMM is available in R (certainly Stata).


text:
How to select point mutations from MAF file
I am using a mutation MAF file from GDC.
I want to select for point mutations, therefore I selected Variant_Classification equal to SNP, and Variant_Classification equal to INS and DEL - but in this case I need to check against the size of the indel mutation.
How can I check for point mutations in the indel case? I think I have to work it out through checking the length of the fields in columns Reference_Allele, Tumor_Seq_Allele1, Tumor_Seq_Allele2 but I am not sure if this is the best solution.


text:
Normalization for two bulk RNA-Seq samples to enable reliable fold-change estimation between genes
I have two bulk RNA-Seq samples, already tpm-normalized.
I would like to know what is a reasonable normalization procedure to enable downstream log fold-change estimation.
The distribution of the two samples using the common set of genes looks similar:

However, the two samples have only been tpm-normalized, is it enough to guarantee reliable fold-change estimation? Should I use another normalization procedure, e.g. Quantile Normalization, before comparison?
My objective is to define a signature using the genes that are up-regulated in Sample1 with respect to Sample0, and vice versa. I'm using log fold-changes, but I'm concerned that their value may be affected by each sample distribution.
Do you also have suggestions for the definition of up-regulated genes with these data?

You have only two samples?
You aren't going to be able to draw strong conclusions from that no matter what you do.  Clever statistics don't work without replicates.

What I have generally done in the past is to process the data using voom in the limma package for bulk RNASeq. Inside voom you can call for different normalization methods to be used - "TMM" works fine for me and, is advocated by many in the field.
voom will output an object containing the normalized expression values in a log2 scale, which, also in my experience, has worked out just fine for calculating log fold changes.
Check out this link for more info on the package as well as normalization methods: https://www.bioconductor.org/packages/devel/workflows/vignettes/RNAseq123/inst/doc/limmaWorkflow.html
It is a very thorough introduction to the package and all of its capabilities.
Good luck!

It's not a good idea to do tpm normalisation prior to differential expression analysis, because the actual read counts are useful to determine shot noise and statistical significance. DESeq2 includes read normalisation as part of its methods for differential expression analysis.
I think shot noise is best explained in terms of shooting photons at a target, as found on Wikipedia:

Shot noise exists because phenomena such as light and electric current
  consist of the movement of discrete (also called "quantized")
  'packets'. Consider light—a stream of discrete photons—coming out of a
  laser pointer and hitting a wall to create a visible spot. The
  fundamental physical processes that govern light emission are such
  that these photons are emitted from the laser at random times; but the
  many billions of photons needed to create a spot are so many that the
  brightness, the number of photons per unit of time, varies only
  infinitesimally with time. However, if the laser brightness is reduced
  until only a handful of photons hit the wall every second, the
  relative fluctuations in number of photons, i.e., brightness, will be
  significant, just as when tossing a coin a few times. These
  fluctuations are shot noise.

Sequencing also has shot noise effects because the sequencing process is a random sampling method. For low-abundance targets, the likelihood of sampling a target is low, so the effect of a single hit is amplified. This reduces the statistical significance when comparing different low-abundance expression values, because the variation in expression is highly dependent on which targets won the sampling lottery. The shot noise variation / expression curve can be seen in Figure 7 of this paper.

I agree with TPM for a lot of reasons, unfortunately the data was already in TPM. Can you explain more about how read counts are useful to determine shot noise and statistical significance? Thanks


text:
Is it possible to use spike-ins with droplet technologies?
By reading Bacher & Kendziorski (2016), I found the following statement:

recent droplet technologies are not yet able to accommodate spike-ins

What's the reason why current droplet technologies cannot use spike-ins? By searching through the available methodologies I found this question on 10X community, but still I cannot understand what is the problem created by using spike-ins with droplet technologies.

Bacher, R. & Kendziorski, C. Design and computational analysis of single-cell RNA-sequencing experiments. Genome Biol. 17, 63 (2016).

I believe it is mostly a cost issue as spike-ins are "expensive" to use given the large number of cells that are/can be sequenced via droplet based methods (as compared to plate-based ones). The 10x link that you provided actually states that this is possible with 10x and Klein et al. used ERCC spike-ins in 953 droplets (Cell. 2015 May 21;161(5):1187-1201. doi: 10.1016/j.cell.2015.04.044.).

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


text:
What's a template switching site?
Reading Islam et al. (2011):

From each transcript, a single read is obtained, corresponding to a template-switching site located preferentially at the 59 end of the mRNA.

By reading this page I think that a template switching site is an additional set of nucleotides added at the 3' end of the newly synthesized cDNA, that serve to anchor a oligo-nucleotide used by the reverse transcriptase to synthesize the second cDNA strand. It is called template switching because reverse transcriptase switches strand upon reaching this set of nucleotides.
Is it correct?

Islam, S. et al. Characterization of the single-cell transcriptional landscape by highly multiplex RNA-seq. genome.cshlp.org
21, 1160–1167 (2011).

From the page you cited:

During first-strand synthesis, upon reaching the 5’ end of the RNA template, the terminal transferase activity of the MMLV reverse transcriptase adds a few additional nucleotide

That usually are CCC (1st is almost always C and the other two most of the times), see: Zajac et al., 2013
These CCCs are used for priming the Template Switching Oligo (TSO), that has GGG and whatever you want to add to the sequence (for instance UMIs and adapters).
So, answering to your question the polimerase does not switch strands, it switches from mRNA as a template to TSO as template. Because of that is called template switching.

You can see it in the graph from the page you cited. In the template switching step polimerase switch from the mRNA to the TSO and continue with the synthesis of the first cDNA strand using TSO as template.


text:
How to get results from Homo.sapiens package in bioconductor for a specific reference
I want to use the Homo.sapiens package in Bioconductor to retrieve the chromosome location start and end for each gene symbol in a specific reference (e.g. hg19 or hg38).
Right now, I am using the following code:
> select(Homo.sapiens, keys=c('KRAS'), cols <-c('SYMBOL', 'TXCHROM', 'TXSTART', 'TXEND'), keytype='SYMBOL')

  SYMBOL TXCHROM  TXSTART    TXEND
1   KRAS   chr12 25358180 25403854
2   KRAS   chr12 25386768 25403863

My question is two-fold:

How can I specify the human genome reference I want to retrieve the data from? (e.g. hg19 or hg38)
Why do I receive multiple and different mappings through Homo.sapiens for some genes? If I check the locations through GeneCards I get:

chr12:25,204,789-25,250,936 (GRCh38/hg38)
chr12:25,357,723-25,403,870 (GRCh37/hg19)

I couldn't find how is this package build, my guess is that is as in the example of the OrganismDBI package. See that vignette section to build your own data package with the genome you want.
The locations are different because you ask for the transcript start and end, if you add the "TXNAME", you will (probably) see there are different transcripts, so different mappings are reported (untested).
Update
The example of the OrganismDBI package builds correctly the Homo.sapiens package, which appears to be linked with hg19. I provide here the relevant code to build the package for hg19 and hg38.
> source("https://bioconductor.org/biocLite.R")
> biocLite("TxDb.Hsapiens.UCSC.hg19.knownGene")
> biocLite("TxDb.Hsapiens.UCSC.hg38.knownGene")

> library(OrganismDbi)

# instructions for hg19, for hg38 change accordingly
> gd <- list(join1 = c(GO.db="GOID", org.Hs.eg.db="GO"),
             join2 = c(org.Hs.eg.db = "ENTREZID",
                       TxDb.Hsapiens.UCSC.hg19.knownGene = "GENEID"))
> destination <- tempfile()
> dir.create(destination)
> makeOrganismPackage(pkgname = "Homo.sapiens.hg19", graphData = gd,
                      organism = "Homo sapiens", version = "1.0.0",
                      maintainer = "Maintainer<maintainer@email>",
                      author = "Author Name", destDir = destination,
                      license = "Artistic-2.0")

Creating package in /var/folders/sf/tfrswpc96jl8b4520mn784440000gn/T//RtmpzR4YnY/file111b490af70f/Homo.sapiens.hg19

# installing package
> install.packages(destination, repos = NULL, type="source")

Then, you can load and use the packages with:
> library(Homo.sapiens.hg19)
> library(Homo.sapiens.hg38)

> select(Homo.sapiens.hg19, keys=c('KRAS'),
     cols = c('SYMBOL', 'TXCHROM', 'TXSTART', 'TXEND'),
     keytype = 'SYMBOL')
'select()' returned 1:many mapping between keys and columns
  SYMBOL TXCHROM  TXSTART    TXEND
1   KRAS   chr12 25358180 25403854
2   KRAS   chr12 25386768 25403863
> select(Homo.sapiens.hg38, keys = c('KRAS'), cols =
       c('SYMBOL', 'TXCHROM', 'TXSTART', 'TXEND'), keytype = 'SYMBOL')
'select()' returned 1:many mapping between keys and columns
  SYMBOL TXCHROM  TXSTART    TXEND
1   KRAS   chr12 25204789 25250931
2   KRAS   chr12 25209168 25250936
3   KRAS   chr12 25209431 25250803
4   KRAS   chr12 25233819 25250929

You're right about TXNAME, I'll check about building my own data package, even if I'd like a tested implementation. Thanks

@gc5 that is the tested implementation, this data packages only collect the information available from other packages. So the data isn't altered

@gc5 Thanks for adding the code you used. I only tweaked the format

Yes, thanks. I still don't understand why the hg19 package does not provide a coherent mapping of transcript start and end with GeneCards, but I think that's minor.


text:
Methodologies for CNV analysis on SNP chip data
I'd like to analyze some Copy Number Variation data produced with SNP chip technologies.
Data is provided using contingency tables, i.e. for each gene it is provided:
         Test_Positive  Test_Negative
Cases           [0,9]+            ...
Controls           ...            ...

In example, test positive may indicate detected heterozygous deletion. In each cell there will be the count of cases/controls tested positive/negative.
Which methodologies may be used other than Fisher's exact test?

There is the hypergeometric test. But what is your end goal, to know if certain case has an association with more CNV or to know the size of the effect of said association?

Is the data format fixed? Because it already implies the test: e.g., having just "cases" and "controls" means that you aren't using any family structures, and I imagine you could have plenty of more ways to group the CNVs rather than just one positive/negative indicator per gene.

@Llopis my end goal is to define the most recurrent and significative copy number variation, intersect them with a set of known loci involved in the disease (which I have in a bedfile) and hopefully run a pathway analysis on the most recurrent or significative ones. However, I don't know if it's sound to use this kind of data for pathway analysis and which tools may be used here

@juod yes the data is unfortunately fixed (at least by now), I have only the indication of single copy amplification/deletion for each gene, and the number of cases/controls affected.

This might be probably answered in stats.stackexchange.com but one of the few alternatives to fisher is the hypergeometric test (but it is closely related to the Fisher's exact test). You could also do a linear model if you tidy up a bit the data. For the pathway analysis I am not sure how to map the CNV and SNP to pathways... but there was a question on the site about it


text:
What is a poly(T) primer anchor sequence?
I am reading Tang et al. (2009).
In the method description it is stated that:

A single cell is manually picked under a microscope and lysed. Then mRNAs are reverse-transcribed into cDNAs using a poly(T) primer with anchor sequence (UP1) and unused primers are digested.

What is an anchor sequence in this case?

Tang, F. et al. mRNA-Seq whole-transcriptome analysis of a single cell. Nat. Methods 6, 377–382 (2009).

As shown in figure 1 of that paper, the anchor sequence is additional sequence that is added to the end of the polyT primer that provides a known and non-target binding site for additional PCR reactions. I gave a presentation about a month ago about the process that Oxford Nanopore Technology use in their strand-switching cDNA kit, which may be useful for understanding the process used by Tang et al.. Here are my strand-switch/PCR images from that presentation:

PolyT primer addition, with mismatched anchor sequence

Strand switch primer addition, with overhanging anchor sequence

PCR amplification, using primers that are complementary to the anchor sequences

In the case of the Tang et al. paper, they use a polyA-tailing method for strand-switching, but the process is essentially the same.

So is it basically a sequence added after the Poly(T) tail, which does not bind to the Poly(A) tail? If so, the main purpose of it is to provide a binding site for the PCR primers?

That's basically correct. It's not "added after", it's part of the same sequence.


text:
Single-cell RNA sequencing (scRNA-seq): filtering cells by transcript counts, how to choose cutoffs?
I am running a notebook with example for the MAGIC algorithm.
In the data preprocessing step, a filtering operation is required to filter out cells with a small count of transcripts.
My question is twofold:

Why filtering out cells with a small count of transcripts?
How to decide the cutoff, i.e. the minimum number of transcripts each cell should have to not be removed from the dataset?

I attach some plots from the original notebook, which display some distributions of the single-cell rnaseq dataset.

In my experience it's helpful to explore the data (ie. make some heatmaps etc) and look for trends that are associated with low molecule counts.  For example, cells with low umi's may look weird or cluster separately.  You can use these to guide your cutoff.

@GWW can you elaborate more on that? Do you mean for example to plot the two principal components and see where the different cells cluster? And maybe if I have a cluster composed by only cells with limited library size, exclude that cluster? Otherwise, which plots will you use? Thanks

I suggest a heatmap of the most variable genes.  PCA may be harder to see.

By reading this tutorial:

Wells with few reads/molecules are likely to have been broken or failed to capture a cell, and should thus be removed.

This answers the first part of the question.
For the second part, i.e. how to decide cutoff, it still seems a bit arbitrary to me, but some authors trim the cells on the left tail of the highest mode, in order to remove other modes (likely due to noise) on the left tail, considering the distribution of cells dependent on the sum of reads (molecules per cell, first plot).


text:
What is the detection limit of protocol in scRNA-seq
In Shapiro et al., when discussing about loss of molecules as source of error in single-cell sequencing, it is written that:

Another source of error is losses, which can be severe. The detection limit of published protocols is $5$–$10$ molecules of mRNA.

I used to think that the detection limit of a protocol was the average maximum number of molecules detected across cells in a scRNA-seq experiment. However, by working on the dataset provided by The Tabula Muris Consortium et al. (pp) I noticed that RNA raw counts often spikes to more than 10 molecules.
What is the actual meaning of detection limit of a protocol in scRNA-seq?

The detection limit refers to the ability to detect a given gene/transcript. In this case, they're saying that a given transcript requires 5-10 copies in order to be likely detected (in other words, the loss rate is 80-90%, which is directly related to a previous question of yours).

Exactly, in that question I misunderstood the concept of detection limit. Now it's clear, thanks.


text:
Which tools are available to call somatic mutations on non-coding regions from whole genome sequencing data?
I have a set of whole genome sequencing samples, some of which have matched normal while some not. I want to call somatic mutations on non-coding regions.
I was looking for GATK best practices (like this) but it seems it is not yet ready for production. Which alternative tools do you suggest?

Have you looked at this question? https://bioinformatics.stackexchange.com/questions/2744/somatic-tumor-only-variant-calling

@Bioathlete thanks, I actually didn't find that question before. However, I found that question too specific on tumor-only samples. I also have normals so maybe there are better tools?

Those will have to be treated as separate analysis.  You can use Mutect or VarScan for the matched samples and the tools listed in the other question for the unmatched samples however that means that you won't be able to compare the results easily.  If the inter-sample comparison is important you will have to run all samples as unmatched.  Personally I have used LoFreq for somatic variant detection with some success.  It is very sensitive so the false positive rate can be high depending on how noisy your data is.


text:
How to compute the chance of failing to detect a gene given the detection limit of a protocol
In Shapiro et al., when discussing about loss of molecules as source of error in single-cell sequencing, it is written that:

Another source of error is losses, which can be severe. The detection limit of published protocols is $5$–$10$ molecules of mRNA. If, as seems likely, the limit of detection is primarily determined by losses during sample preparation, this would indicate that $80$–$90\%$ of mRNA was lost. Or, to put it the other way around, a $90\%$ loss leads to an approximately $50\%$ chance of failing to detect a gene that is expressed at a level of seven mRNA molecules (from the binomial distribution).

How is this probability computed using the binomial distribution? I thought that $90\%$ loss corresponds to $5$ detected molecules, and I assume that $k=7$ for the binomial calculation, but I am unable to go further.

A 90% loss can be rephrased as a 10% chance of detecting anything. So what we want to find is the probability of detecting 0 molecules, when we start with 7 and have 10% probability of success. Once can do that in R as follows:
> pbinom(0, 7, 0.1)
0.4782969

So ~50%, as they stated. I suspect that part of the confusion arises from the fact that the detection limit is due to lowly-expressed genes/transcript being heavily affected by this loss. So the probability of detecting a single molecule out of 4 original molecules (assuming 90% mRNA loss) is 34%, for 3 molecules it's 27%, for 2 it's 19% and for 1 it's 10%. I think the threshold of 5 is mentioned more because it's a nice round number than there's anything particularly different in detecting a gene with 4 vs. 5 molecules (34 vs. 41% probability).

I was reading again your answer to understand the first sentence. If I understood correctly, the mapping from $90%$ loss to $10%$ chance happens computing the probability of detecting a single molecule out of $1$ original molecule. Right?

Correct, since there's a 90% loss of all signal, there's a 10% probability of detecting any given molecule.

Detecting zero molecules out of 7 can be rephrased as failing to detect a molecule, 7 times independently. The probability of this is 0.9 to the power of 7. Python tells me (0.9 ** 7) that this is 0.4782969. That's probably what the modelization in terms of binomial distribution means, but to me it is much easier to understand in terms of this very basic probability calculation.


text:
What is local realignment and what is the problem it solves?
I am trying to understand local realignment but I could not get a clear idea of what is the problem solved by it.
For example, reading Homer and Nelson (2010):

Because alignment algorithms map each read independently to the reference genome, alignment artifacts could result, such that SNPs, insertions, and deletions are improperly placed relative to their true location. This leads to local alignment errors due to a combination of sequencing error, equivalent positions of the variant being equally likely, and adjacent variants or nearby errors driving misalignment of the local sequence. These local misalignments lead to false positive variant detection, especially at apparent heterozygous positions. For example, insertions and deletions towards the ends of reads are difficult to anchor and resolve without the use of multiple reads. In some cases, strict quality and filtering thresholds are used to overcome the false detection of variants, at the cost of reducing power [13]. Since each read represents an independent observation of only one of two possible haplotypes (assuming a diploid genome), multiple read observations could significantly reduce false-positive detection of variants.

I don't understand how the different artifacts (such as SNPs) can result from the alignment.

Homer, Nils, and Stanley F. Nelson. "Improved variant discovery through local re-alignment of short-read next-generation sequencing data using SRMA." Genome biology 11.10 (2010): R99.

SNPs are likely to be created and InDels are likely to be missed. Suppose you have a read, ACTGACTGACTGTAC and you align it to a reference sequence ACTGACTGACTGTTAAGAACGACTACGAC. If you aligned that, you would either get:
ACTGACTGACTGTac    (lower case denotes soft-clipping)
ACTGACTGACTGTTAAGAACGACTACGAC

or
ACTGACTGACTGTAC    (N.B., you've created some possible SNPS)
ACTGACTGACTGTTAAGAACGACTACGAC

But what happens if you look at a bunch of reads? See below:
ACTGACTGACTGTAC
  TGACTGACTGTACTA
    ACTGACTGTACTAAG
ACTGACTGACTGT--TAAGAACGACTACGAC

I've made the 2 base insertion obvious. As you can see, using many alignments, we can "realign" our original read to see that it actually ends in a 2 base insertion, something that's completely impossible to find with a single read. This realignment is part of what the GATK haplotype caller does (it makes a local de Bruijn graph). That's also what my MethIndelRealigner does. I have another real-world example in its documentation.
Raw alignments in IGV:

Alignments after realignment:

One can find plenty of similar scenarios that result in SNPs rather than simply losing coverage for InDels.

Awesome.. I was not considering the problem of mapping each read independently. Now it's clear, thanks!


text:
taxon exclude list for searching local blast database using blastn
I am looking for a solution to exlude certain entries when searching a local blast nt database (with blastn), specifically the sequences from uncultured / environmental samples, ideally using their taxon ids.
This should ideally replicate the results from the BLAST web interface, when checking the option "Exclude Uncultured/environmental sample sequences", see screenshot below.
Three approaches that I am aware of so far are:
(1) grep taxon IDs for any taxa with "uncultured" or "environmental" in their scientific name from names.dmp and fetch the GIs belonging to these taxa using the E-utils. Then use the -negative_gilist parameter for blastn. However, GIs are discontinued for future records, so this method should be replaced by Acc.Ver identifiers.
(2) locally format the blast database after prepocessing the nt file by myself and exclude records belonging to the taxa list as obtained in (1).
(3) Postprocess the blast output and remove the relevant lines, but this approach might remove a lot of alignments, whereas (1) and (2) would always generate the desired number of output alignments.
So, I am looking for other suggestions, preferably without preprocessing the nt file.
Ideally, one would like to have a -negative_taxalist option for blastn, which would resemble what is possible on the NCBI BLAST web interface, where one can exclude taxa, see screenshot below.

Thanks!

Hurray, the solution seems already to be implemented in the new BLAST database format (v5), which is usable with blast+ versions >= 2.8.0.
It is possible to specify an exclusion as well as inclusion list of taxon ids, which, however, have to be on species level.
Species level taxon ids can be obtained from higher level taxon ids using the script get_species_taxids.sh, e.g.: get_species_taxids.sh -n Enterobacterales
See docs here: https://ftp.ncbi.nlm.nih.gov/blast/db/v5/blastdbv5.pdf
Edit:
For creating taxon ID lists locally from the NCBI taxonomy files, the program subtree can be used. See example therein for compiling a list of taxon IDs from environmental or unclassified taxa.

Sorry to bother you, but how did you manage to work with get_species_taxids.sh? Because I installed edirect on linux, but won't work with the script.


text:
What is "unmapped read segments" in the output of samtools idxstats?
samtools idxstats produces a four column output  (see here)
The output is TAB-delimited with each line consisting of reference sequence name, sequence length,
# mapped read-segments and # unmapped read-segments. It is written to stdout. Note this may count reads
multiple times if they are mapped more than once or in multiple fragments.

What is the meaning of unmapped reads here? Naively, I would take it to mean all the reads that didn't map to a specific feature... however in a sample with millions of reads most of them should be counted as such, which is not the case in my calculations, where the number of unmapped reads is of the same order as the number of mapped reads.
I will very much appreciate explanations from people experience with samtools.

An unmapped read segment is any entry with bit 4 of the flag field set (i.e., it's unmapped). Assuming the reference to which you aligned reasonably matches the sample from which the reads were generated, less than ~10% of to total reads should be unmapped.
Please note that BAM files have no idea what a "feature" is. They only contain information about reads and the chromosomes/contigs to which they align.

I have mapped to a gene catalogue, i.e., I have many references, and every line in the idxstats output correspondes to one such reference. I understand how a read mapped to a reference/gene/feature is associated with it, but a read unmapped is unmapped to all of them...
But what you wrote about flags makes a lot of sense - I have to look what happened at the mapping step. Thanks.

For paired-end reads, some aligners will assign an unmapped read to the chromosome of its mate if its mate aligned (the unmapped bit in the flag will still be set).


text:
Reads count in metagenomics
Background:  I am developing a pipeline for metagenomic studies of human gut microbiote. In particular, I am mapping the reads data originating from shotgun whole genome sequencing onto a gene catalogue (similar to IGC) and counting reads mapped on each feature/gene. Till now I have been using NGLess, which is a rather convenient tool. However, we have some questions about how exactly it counts reads and it might be not the best option for scaling up our pipeline.
Question:
I am looking for alternative tools for gene count. A brief search has shown that htseq-count and featureCount are two frequently used tools. However they mainly appear in RNAseq context, and I am not certain whether they could be used in my case. I would also appreciate comments on the possibility of using samtools mpileup.

I think that there is a good argument for using technologies adapted from RNA-seq for metagenomic abundance quantification. The tool that I am most familiar with (that I use for this purpose) is kallisto, which uses pseudoalignment instead and is therefore less compute intensive.
For some background on the conceptual linkages between metagenomic abundance and transcriptomic abundance you can see the kallisto metagenome abundance paper or the explanatory blog post.
I think that it is preferable to use a tool that is specifically adapted for metagenomics, as there are always some numerical issues that might have been accounted for by the authors. I think that for this purpose it is maybe best not to naively use tools like mpileup, though I have no specific objection to it. However, if you look around you can see various papers using e.g. Cufflinks and Cuffdiff for effectively the same purpose, and Cufflinks is definitely intended for RNA-seq.
Depending on your specific application, you might want to use a different tool like mOTUS2, if you are interested in estimating abundance of specific genomes or lineages in a metagenome mixture. This method is probably more accurate as it is based on averaging across the whole genome, but it is not very useful if what you are interested in is more e.g. functional profiling of the metagenome.

Thank you for answering. We did try motus2, but it is restricted to a custom-made gene database and the marker genes grounded in the associeted species classification. What I am looking for is a tool for creating a gene count table, which we then could use to develop our own methods of species identification.

@Vadim In that case, I think that a tool like kallisto (or even cufflinks) might be most helpful. These will be providing counts as TPMs (though they also provide raw counts), of course, so if there is some statistical issue with that kind of compositional count then it might be a problem. However, metagenomics is intrinsically compositional so I'm not sure other methods will be much better. If you are ultimately going to do deconvolution of some kind, it might be helpful to look at what metabat2 is doing (more for ideas than for usefulness as a tool, as it's solving a different problem).


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

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

This is however complicated by such intricacies of sequencing as:

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

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

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

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


text:
Comparing two genome annotations
I have one-dimentional array (human genome). Also I have two annotations for it, we can think about them as different peaks (it's nucleosome and secondary structures). How can we find correlation and/or causation between this annotations? Currently we are forming hypothesis like "one element is always before other" and statistically testing them. Can we be smarter and use known tools or at least automate this?

Do you mean with something like a Markov model?

Yeah, but how will you implement them?

I am also interested in applying such models to multiple 'omics data, I hope this question will get some useful answers.

Could you specify which kind of hypothesis do you want to test? What I have seen is centering all the positions you are interested in and observe if there is a common pattern on the other variables observed (annotations, CpG, mutation rates...) which is not noise. See here an example

Perhaps phrase this as time series analysis? Does one event precede another in a statistical sense? I think if you look for tools in this domain, you'll find some good solutions.

See answers to this plot https://bioinformatics.stackexchange.com/questions/83/how-to-correlate-two-zero-inflated-bedgraph-like-signals?rq=1

I would start with metaplots for a visual/exploratory analysis to guide your statistical testing.  Cut the genome into short segments e.g. centred on the centre of each nucleosome peak, overlay all the segments on top of each other, and plot the abundance of secondary structure features at each base count as you step away from the centres.  This should give an easy view on whether there is any clear structure in the relative positions of the different objects.
Do the plots both ways round, and try different subsets of the objects (e.g. ones overlapping genes or not, different stringency thresholds for peak calling).  This should hopefully give you some insights on what is worth testing for.

You can express the annotations of the genome with a dot plot, x and y represent the two annotations along the genome. The you can summary the result with Spearman correlation coefficient(s) based on whole set of annotations or particular gene families.


text:
What are Approximate Read Counts (Library Sizes) and Lengths (Insert Sizes) for Next-Generation DNA Sequencers?
I am performing a simulation study and am curious about the parameters of my simulated metagenome.
What are the library and insert sizes of some of the most "used" sequencers. Mainly, I am interested in the different Illumina platforms (MiSeq, HiSeq, HiSeqX), 454 Pyrosequencing, and Ion Torrent.
I do not need the long read technologies like Pac-Bio's SMRT and Oxford Nanopore since they are not typically used with metagenomics, yet.

The library size depends usually on your budget. There are some reviews where they compare the different platforms, but these papers get outdated so quickly, that it would be better to just check wiki for the specs of different platforms. Also on the websites of Illumina or Life Tech, you'll find the specs.
NB. 454 sequencing is discontinued, so that is not very relevant anymore in my opinion.

Thanks for the Wiki link. This is what I was looking for!

Both of those depend on the application of the library.  For example whole genome assembly of a mammal would require a larger library size than targeted re-sequencing of a virus.  Also insert size varies if you are doing genome assembly a larger insert size is valuable to allow for long range linking and scaffolding.
I think that if you edit your question of be more specific people will be able to provided a more concrete answer other wise it is "it depends"


text:
Converting hmmer --tblout to pandas dataframe
Is there a way to convert a hmmer output to a pandas dataframe? I am also unsure how to load a hmmer tblout table into python via the Bio module. I believe you can call a hmmer format with SeqIO.parse or SeqIO.search. How can you convert the table to a pandas dataframe?

This question was answered on Stack Overflow: https://stackoverflow.com/questions/62012615/convert-a-hmmer-tblout-output-to-a-pandas-dataframe/62021471#62021471

It is rude to cross-post on multiple forums.

The manual says fields in a tabular output file are space-delimited. It should be fairly straightforward to load a space-delimited file in pandas.
From the manual:

--tblout <f> Save a simple tabular (space-delimited) file summarizing the per-target output, with one data line per homologous target model found.

While this link may answer the question, it is better to include the essential parts of the answer here and provide the link for reference.  Link-only answers can become invalid if the linked page changes. - From Review

Thank you for the pointer, @Michael. I've added an excerpt from the manual.

Thanks @RamRS ... much better :-)

This question was answered on StackOverflow:  https://stackoverflow.com/questions/62012615/convert-a-hmmer-tblout-output-to-a-pandas-dataframe/62021471#62021471


text:
Transferring genomic features on new coordinates
I have a eukaryotic genome for which an updated sequence for a chromosome  was recently obtained. I want to map RNAseq reads on the genome (and perform other downstream analyses) and would like to use the most up-to-date information possible (so the 'new' sequence of that chromosome).
However, before performing the read mapping, I'd like to update the annotation of my genome (GFF3) to match with the 'new' coordinates of this chromosome.
Basically, I would like to transfer information from the 'old-coordinates' GFF3 to obtain a 'new-coordinates' GFF3, and if possible conserve all the information/hierarchy of the file (gene, mRNA, exon, etc.).
Example:
chrXX   source  gene    222 5942    .   -   .   ID=gene_1;Name=gene_1;length=5720
chrXX   source  mRNA    222 5942    .   -   .   ID=gene_1.1;Parent=gene_1;Name=gene_1.1;length=5720
chrXX   source  exon    222 5794    .   -   .   ID=gene_1.1.2;Parent=gene_1.1
chrXX   source  exon    5889    5942    .   -   .   ID=gene_1.1.1;Parent=gene_1.1
chrXX   source  CDS 222 5794    .   -   1   ID=CDS:gene_1.1.2;Parent=gene_1.1;Name=gene_1.1
chrXX   source  CDS 5889    5942    .   -   0   ID=CDS:gene_1.1.1;Parent=gene_1.1;Name=gene_1.1

... should be updated to ...
chrXX   source  gene    333 6053    .   -   .   ID=gene_1;Name=gene_1;length=5720
chrXX   source  mRNA    333 6053    .   -   .   ID=gene_1.1;Parent=gene_1;Name=gene_1.1;length=5720
chrXX   source  exon    333 5905    .   -   .   ID=gene_1.1.2;Parent=gene_1.1
chrXX   source  exon    6000    6053    .   -   .   ID=gene_1.1.1;Parent=gene_1.1
chrXX   source  CDS 333 5905    .   -   1   ID=CDS:gene_1.1.2;Parent=gene_1.1;Name=gene_1.1
chrXX   source  CDS 6000    6053    .   -   0   ID=CDS:gene_1.1.1;Parent=gene_1.1;Name=gene_1.1

The approach I tried (but I am not sure it is the way to do it):

Extract sequence of all the features that are going to be updated to fasta (bedtools getfasta)
Map these sequences to the genome with the updated chromosome (gmap with --nosplicing since sequences we are mapping correspond to genomic regions).
Create the updated GFF3 file. When doing that, I apply one supplementary rule: if a feature mapped on another chromosome, but there was also another alignment on the updated chromosome, prioritize the one that correspond to the same (updated) chromosome.

What would be the proper way to do such a thing?  Any suggestion of methods/tools is welcome!

I'd suggest to paste the full code used for each step, so people can advise you better.

I think the standard way of doing this is to make a chain file, then use it to do a liftOver of the annotations:
mkdir psl
for i in ../ci3/rm/masked/*.masked; do blat ../ci2.2bit $i -tileSize=12 -fastMap -minIdentity=98 psl/`basename $i .fa.masked`.psl -noHead -minScore=100; done

Translate psl files to chains in the directory chain:
mkdir chain
for i in psl/*.psl; do axtChain -linearGap=medium -psl $i ../ci2.2bit ../ci3/ci3.2bit chain/`basename $i .psl`.chain; done

Merge short chains into longer ones into the directory chainMerge:
mkdir chainMerge
chainMergeSort chain/*.chain | chainSplit chainMerge stdin -lump=50

concat and sort the chains:
cat chainMerge/*.chain > all.chain
chainSort all.chain all.sorted.chain

Need info about chromosome sizes for netting:
twoBitInfo ../ci3/ci3.2bit ci3.chromInfo
twoBitInfo ../ci2.2bit ci2.chromInfo

Netting: identify alignable regions from chains:
mkdir net
chainNet all.sorted.chain ci2.chromInfo ci3.chromInfo net/all.net /dev/null

Finally, select the right alignable regions using the nets, creating a "liftOver" file:
netChainSubset net/all.net all.chain ci2ToCi3.liftOver

Run liftOver:
CrossMapy.py bed ci2ToCi3.liftOver test.hg18.bed

Thanks for the quick answer, will test it and report back.

In addition to the answer by heathobrien you might take a look at RATT: Rapid Annotation Transfer Tool. The approach is essentially the same, only using more modern approaches to for instance whole genome alignments. It uses nucmer, which when using the latest version (4.0) allows you to do whole genome alignments using parallel computing.

Are there any resources that show how to do this? The documentation seems really non-user friendly.


text:
What do 0s in stageR mean?
I'm running stageR as part of a pipeline based on DRIMSeq and wanted to know what it means when the transcript level padj is equal to 0, particularly when the gene is non-significant. I'm pretty sure that in some cases if the padj is 0 it means that the chance of it being different due to chance is 0. There are a lot more 0s than I would expect in sites where the gene padj is non-significant. One example, I have a gene padj that's >0.9 but both of the sites I'm looking at have a padj of 0.
I have a couple of suspicions about why this is: either because it defaults to this if the comparative difference is or it uses 0 as a substitute to NA. Am I thinking along the right lines? Are the 0s real or are they likely to be a quirk of the package?

In stageR, we automatically return stage-wise adjusted p-values of zero for transcripts coming from genes that were found to be significant and only have 1 transcript, when testing for differential transcript expression (DTE).
Similarly, when testing for differential transcript usage (DTU), we return p-values of zero for transcripts derived from a significant gene consisting of no more than two transcripts (based on the tx2gene object you provide).
Transcript-level p-values for non-significant genes should not occur. We provide NA's for genes with only one transcript for DTU, and by default you should get NA p-values for non-significant genes, since these genes do not proceed further than the screening stage.
Could you provide some more details on the analysis? Which multiple testing correction have you selected?
Would it be possible to provide a reproducible example for your analysis?

I'm using DTU and I've set the alpha to 0.5 as we were hoping that some of the genes which were non-significant overall might have some features that were significant. I'll see if I can make a reproducible example in a bit.

I'm really sorry it took so long to get back to you, I ended up getting distracted by something else. It turns out that the issue was probably elsewhere within what I was doing.

No problem. Let us know if you have any remaining questions or issues.


text:
Best way to align to find inserted sequence
We have some RNA from knock-in mice, there are two different sequences we're looking for. We have aligned to the mouse genome using STAR but the sequence isn't there which isn't too surprising What is the best way to check for the presence of the sequence? Should I just turn them into a sequence and use them as the reference genome for STAR to align to? Is there a way I can easily incorporate it into the existing genome so I could get read counts for everyhting including those reads? Should I just BLAST using something like magicblast?

Depends on the throughput. High throughput download the mouse genome and use Blat (a bit out of date now)? Megablast just make sure you select the mouse genome to speed the blasting.

@ATpoint could you give me some suggestions on what to search by any change? I couldn't really find anything useful partly because it kept talking about Indels.

@Michael I'll look into that. Thanks.


text:
Cross-reference with PDB database
I have a list of several thousand proteins and their UNIPROT IDs. I'm looking for an efficient method of cross-referencing it against the PDB tertiary structure database, and get a list of those proteins with a tertiary structure in the PDB database.
I've tried to BLASTP the list of UNIPROT IDs against the PDB database, using the NCBI BLAST portal but encountered too many errors of "Error: Failed to read the Blast query: Sequence ID not found", making the process of manual filtering not convenient and not efficient.
Any ideas?
Thank you!

You can use the SIFTS mapping: http://www.ebi.ac.uk/pdbe/docs/sifts/
From an article:

The Structure Integration with Function, Taxonomy and Sequences
  resource (SIFTS; http://pdbe.org/sifts) is a close collaboration
  between the Protein Data Bank in Europe (PDBe) and UniProt. The two
  teams have developed a semi-automated process for maintaining
  up-to-date cross-reference information to UniProt entries, for all
  protein chains in the PDB entries present in the UniProt database.
  This process is carried out for every weekly PDB release and the
  information is stored in the SIFTS database.

It has no web interface for querying. You just download CSV (or TSV) file with all the mapping data.

Thank you. I'll try it if the ID mapping tool won't provide.

The UniProt website offers an ID mapping tool that tends to work pretty well:

Then, hit "submit" and you will get the list of corresponding PDB IDs:

Thank you, it worked.
Please aid a PDB-Noob like myself to understand the results: https://ibb.co/jTMGmb
How come each of the UNIPROT ID received several different PDB IDs?
https://ibb.co/b5fRmb

@AdrianSmith that's how PDB works. You'll have one structure for the active site, another for the active site bound by foo, one for the active site bound by bar, another structure for the protein in complex with something else. . . This should be clear enough if you just click on the results you got and see their descriptions.

Thanks for the explanation. I'm looking for a specific protein-target structure (I have a list of ligand-protein interactions and i'm looking for the 3d structure of the binding site if available). Is there any way to search through all of the different "structures"? Does the PDB site offer such a curated database?

@AdrianSmith no idea, I'm afraid. I've almost never used PDB. I suggest you ask a new question explaining what you need in more detail. Ideally, give an example UniProt ID and explain what you would want to see from it.


text:
Using R to import specific records from a PDB file
I'm using R and I'm trying to import specific records from a PDB file (e.g. HETNAM, HETSYN). Unfortunately, it seems the read.pdb function of either biomaRt or Rpdb aren't capable of it.
Any ideas?
Thanks

Hi, could you show what have you tried? As far as I know it is not so difficult to parse the files, but that depends on the format of the file you downloaded.

Hi, nothing special, I just used the read.pdb function to read in the PDB files.

If your file is in PDB format as the one shown here in .pdb format then you can simply extract your desired fields using sed which is a programming language used to parse and transform text. As for this structure i.e 6C4A in .pdb format HETNAM and HETSYN can be extracted using sed as follows:
sed  -n -e '/HETNAM/ p' -e '/HETSYN/ p'  6C4A.pdb

HETNAM     EJA S-[(1Z)-2-CARBOXY-N-HYDROXYETHANIMIDOYL]-L-CYSTEINE
HETNAM     PYR PYRUVIC ACID
HETNAM      MG MAGNESIUM ION
HETNAM     3NP 3-NITROPROPANOIC ACID
HETNAM     PEG DI(HYDROXYETHYL)ETHER
HETNAM     EPE 4-(2-HYDROXYETHYL)-1-PIPERAZINE ETHANESULFONIC ACID
HETSYN     EPE HEPES

Thanks, seems pretty straightforward. How do I use sed? Is it a way to apply it on R ?

It is straight forward. Write this command on Linux command line/terminal with your input file present in current working directory.

Secondly for using sed in R you should check this


text:
How can I create my own GO slim?
I know there are existing GO slims predefined. But does anyone know how I can create my own collection of GO terms for enrichment analysis?

I have never created one, but a GO slims is just a list of GO associated with genes, so you just need to download the list of GO and associate them with your list of genes.  What have you tried/searched? What do you need that it is not covered in the FAQ?

GO slims are just a subset of the gene ontology graph. Just a group of terms, no genes. We can't help you create one if you don't explain what terms you want and what sort of analysis you want to do.

I have a list of genes enriched in kidney. So I want to get a list of kidney-related terms and see which of them are significantly enriched.

@Dr_Hope Then what you really want is to test an enrichment in GO related to kidney  not create a GO slim. But if you only use a subset of genes and a subset of GOs the test would be biased (AFAIk) if you don't tell the suitable background population of genes and GO terms.

I've already done a canonical GO to check the background. The enriched terms are mostly kidney-related because our sample is a kidney. That's why I need a GO slim to check kidney-only terms.

authors: กิ๊ฟเลิกล่ะสายแข็ง


text:
How can I download genePred file from UCSC?
I need a gtf annotation file. UCSC doesn't give us a proper gtf file with distinct gene_id and transcript_id. It asked us to get a genePred file to convert to gtf.
But where can we get genePred files?

It turned out that using table browser on UCSC we could download an output with all fields. That output actually contains a GenePred file but has one extra column on the left end.
The genePred file format definition is here.


text:
CIBERSORT runtime error
I was running CIBERSORT web and saw this error:
Exception in thread "Thread-0" java.lang.ArrayIndexOutOfBoundsException: 0
at launch.QN.doQN(QN.java:56)
at launch.Cibersort.Run_CIBERSORT(Cibersort.java:611)
at launch.Cibersort.run_all(Cibersort.java:179)
at launch.Cibersort$runCIBERSORT.run(Cibersort.java:1585)
at java.lang.Thread.run(Thread.java:745)
311

Does anyone know why?

OK I found the bug. CIBERSORT seems to require tsv files instead of csv files. Plus it doesn't seem to like empty columns