Genome: Difference between revisions
Jump to navigation
Jump to search
| Line 787: | Line 787: | ||
=== curatedTCGAData === | === curatedTCGAData === | ||
* [https://www.bioconductor.org/packages/release/data/experiment/html/curatedTCGAData.html Curated Data From The Cancer Genome Atlas (TCGA) as MultiAssayExperiment Objects] | |||
** data access through ExperimentHub | |||
** provides data from the legacy Firehose pipeline | |||
* [https://waldronlab.io/MultiAssayWorkshop/ Multi-omic Integration and Analysis of cBioPortal and TCGA data with MultiAssayExperiment] from [https://bioc2020.bioconductor.org/workshops.html Bioc2020] | * [https://waldronlab.io/MultiAssayWorkshop/ Multi-omic Integration and Analysis of cBioPortal and TCGA data with MultiAssayExperiment] from [https://bioc2020.bioconductor.org/workshops.html Bioc2020] | ||
* [https://waldronlab.io/PublicDataResources/ Public data resources and Bioconductor] from [https://bioc2020.bioconductor.org/workshops.html Bioc2020] | * [https://waldronlab.io/PublicDataResources/ Public data resources and Bioconductor] from [https://bioc2020.bioconductor.org/workshops.html Bioc2020] | ||
Revision as of 13:27, 24 August 2020
Visualization
See also Bioconductor > BiocViews > Visualization. Search 'genom' as the keyword.
IGV
- UserGuide
- Inspect alignments with IGV in Anders2013
- For paired data we will see two reads (one is -> and the other is <- direction) from the same fragment.
nano ~/binary/IGV_2.3.52/igv.sh # Change -Xmx2000m to -Xmx4000m in order to increase the memory to 4GB ~/binary/IGV_2.3.52/igv.sh
Simulated DNA-Seq
The following shows 3 simulated DNA-Seq data; the top has 8 insertions (purple '|') per read, the middle has 8 deletions (black '-') per read and the bottom has 8 snps per read.
Whole genome
Whole exome
- (Left) GSE48215, UCSC hg19. It is seen there is a good coverage on all exons.
- (Right) 1 of 3 whole exome data from SRP066363, UCSC hg19.
File:Igv gse48215.png File:Igv srp066363.png
RNA-Seq
- (Left) Anders2013, Drosophila_melanogaster/Ensembl/BDGP5. It is seen there are no coverages on some exons.
- (Right) GSE46876, UCSC/hg19.
File:Igv anders2013 rna.png File:Igv gse46876 rna.png
Tell DNA or RNA
- DNA: no matter it is whole genome or whole exome, the coverage is more even. For whole exome, there is no splicing.
- RNA: focusing on expression so the coverage changes a lot. The base name still A,C,G,T (not A,C,G,U).
RNA-seq DRaMA
https://hssgenomics.shinyapps.io/RNAseq_DRaMA/ from 2nd Annual Shiny Contest
Gviz
GIVE: Genomic Interactive Visualization Engine
Build your own genome browser
ChromHeatMap
Heat map plotting by genome coordinate.
ggbio
NOISeq package
Exploratory analysis (Sequencing depth, GC content bias, RNA composition) and differential expression for RNA-seq data.
rtracklayer
R interface to genome browsers and their annotation tracks
- Retrieve annotation from GTF file and parse the file to a GRanges instance. See the 'Counting reads with summarizeOverlaps' vignette from GenomicAlignments package.
ssviz
A small RNA-seq visualizer and analysis toolkit. It includes a function to draw bar plot of counts per million in tag length with two datasets (control and treatment).
Sushi
See fig on p22 of Sushi vignette where genes with different strands are shown with different directions when plotGenes() was used. plotGenes() can be used to plot gene structures that are stored in bed format.
cBioPortal and TCGA
- The cBioPortal for Cancer Genomics provides visualization, analysis and download of large-scale cancer genomics data sets.
- https://github.com/cBioPortal/cbioportal
- https://gdac.broadinstitute.org/
- TCGA processed RNA-Seq data (GSE62944) as a SummarizedExperiment
- Tutorial: retrieve full TCGA datasets from cBioportal with R
- cgdsr package
- 生物信息学生R入门教程
Qualimap
Qualimap 2 is a platform-independent application written in Java and R that provides both a Graphical User Inteface (GUI) and a command-line interface to facilitate the quality control of alignment sequencing data and its derivatives like feature counts.
SeqMonk
SeqMonk is a program to enable the visualisation and analysis of mapped sequence data.
Copy Number
Copy number work flow using Bioconductor
Detect copy number variation (CNV) from the whole exome sequencing
- http://www.ncbi.nlm.nih.gov/pubmed/24172663
- http://www.biomedcentral.com/1471-2164/15/661
- http://www.blog.biodiscovery.com/2014/06/16/comparison-of-cnv-detection-from-whole-exome-sequencing-wes-as-compared-with-snp-microarray/
- Bioconductor - CopyNumberVariation
- exomeCopy package
- exomeCNV package and more info including slides by Fah Sathirapongsasuti.
- ximmer and source
- dudeML: A Simple Deep Learning Approach for Detecting Duplications and Deletions in Next-Generation Sequencing Data
Whole exome sequencing != whole genome sequencing
Consensus CDS/CCDS
- UCSC
- http://genome.ucsc.edu/cgi-bin/hgTrackUi?hgsid=573539327_9Ga3TjKR7LNFYJ2JIMZKDNZfl9eF&c=chr1&g=ccdsGene
- Schema/Example
- The GenomicFeatures package can download the CCDS regions.
- The CCDS regions are convenient to use as genomic ranges for CNV detection in exome sequencing data, as they are often in the center of the targeted region in exome enrichment protocols and tend to have less variable coverage then the flanking regions -- exomeCopy package vignette.
DBS segmentation algorithm
DBS: a fast and informative segmentation algorithm for DNA copy number analysis
modSaRa2
An accurate and powerful method for copy number variation detection
NGS
File:CentralDogmaMolecular.png
See NGS.
R and Bioconductor packages
Resources
- HarvardX Biomedical Data Science Open Online Training, Bioinformatics Training at the Harvard Chan Bioinformatics Core
- CSAMA 2017: Statistical Data Analysis for Genome-Scale Biology
- Some basics of biomaRt (and GenomicRanges)
- Annotating Ranges Represent common sequence data types (e.g., from BAM, gff, bed, and wig files) as genomic ranges for simple and advanced range-based queries.
library(VariantAnnotation) library(AnnotationHub) library(TxDb.Hsapiens.UCSC.hg19.knownGene) library(TxDb.Mmusculus.UCSC.mm10.ensGene) library(org.Hs.eg.db) library(org.Mm.eg.db) library(BSgenome.Hsapiens.UCSC.hg19)
- Analysis of RNA-Seq Data with R/Bioconductor by Girke in UC Riverside
- 2018 Data Intensive Biology Summer Institute at UC Davis
- R / Bioconductor for High-Throughput Sequence Analysis 2012 by Martin Morgan1 and Nicolas Delhomme
- Count-based differential expression analysis of RNA sequencing data using R and Bioconductor by Simon Anders
- Sequences, Genomes, and Genes in R / Bioconductor by Martin Morgan 2013.
- Orchestrating high-throughput genomic analysis with Bioconductor by Wolfgang Huber et al 2015.
- RNA-seq Analysis Example This is a script that will do differential gene expression (DGE) analysis for RNA-seq experiments using the bioconductor package edgeR. RPKMs were calculated for bar plots.
- An introduction to R and Bioconductor for the analysis of high-throughput sequencing data by Pascal MARTIN Oct 2018
Docker
Bioinstaller: A comprehensive R package to construct interactive and reproducible biological data analysis applications based on the R platform. Package on CRAN.
Some workflows
RNA-Seq workflow
Gene-level exploratory analysis and differential expression. A non stranded-specific and paired-end rna-seq experiment was used for the tutorial.
STAR Samtools Rsamtools
fastq -----> sam ----------> bam ----------> bamfiles -|
\ GenomicAlignments DESeq2
--------------------> se --------> dds
GenomicFeatures GenomicFeatures / (SummarizedExperiment) (DESeqDataSet)
gtf ----------------> txdb ---------------> genes -----|
RNA-Seq differential expression work flow using DESeq2
Sequence analysis
library(ShortRead) or library(Biostrings) (QA) gtf + library(GenomicFeatures) or directly library(TxDb.Scerevisiae.UCSC.sacCer2.sgdGene) (gene information) GenomicRanges::summarizeOverlaps or GenomicRanges::countOverlaps(count) edgeR or DESeq2 (gene expression analysis) library(org.Sc.sgd.db) or library(biomaRt)
Accessing Annotation Data
Use microarray probe, gene, pathway, gene ontology, homology and other annotations. Access GO, KEGG, NCBI, Biomart, UCSC, vendor, and other sources.
library(org.Hs.eg.db) # Sample OrgDb Workflow
library("hgu95av2.db") # Sample ChipDb Workflow
library(TxDb.Hsapiens.UCSC.hg19.knownGene) # Sample TxDb Workflow
library(Homo.sapiens) # Sample OrganismDb Workflow
library(AnnotationHub) # Sample AnnotationHub Workflow
library("biomaRt") # Using biomaRt
library(BSgenome.Hsapiens.UCSC.hg19) # BSgenome packages
| Object type | example package name | contents |
|---|---|---|
| OrgDb | org.Hs.eg.db | gene based information for Homo sapiens |
| TxDb | TxDb.Hsapiens.UCSC.hg19.knownGene | transcriptome ranges for Homo sapiens |
| OrganismDb | Homo.sapiens | composite information for Homo sapiens |
| BSgenome | BSgenome.Hsapiens.UCSC.hg19 | genome sequence for Homo sapiens |
| refGenome |
RNA-Seq Data Analysis using R/Bioconductor
- https://github.com/datacarpentry/rnaseq-data-analysis by Stephen Turner.
- Tutorial: Introduction to Bioconductor for high-throughput sequence analysis by UseR 2015
- systemPipeR Building end-to-end analysis pipelines with automated report generation for next generation sequence (NGS) applications such as RNA-Seq, ChIP-Seq, VAR-Seq and Ribo-Seq. An important feature is support for running command-line software, such as NGS aligners, on both single machines or compute clusters.
- BioC2015
recount2
- https://github.com/leekgroup/recount
- recount2 - A multi-experiment resource of analysis-ready RNA-seq gene and exon count datasets.
- Human RNA-seq data from recount2 and related packages from BioC 2020 Workshop
- recount2 works with DESeq2, edgeR and limma-voom
dbGap
eQTL
Statistics for Genomic Data Science (Coursera) and MatrixEQTL from CRAN
GenomicDataCommons package
- Genomic Data Commons
- Use the GenomicDataCommons package to find and download variants from TCGA (NCI Genomic Data Commons Access) dataset and maftools package for analysis and visualization. See https://seandavi.github.io/talk/2018/02/08/bioconductor-a-potential-hub-in-the-cancer-biomarker-data-ecosystem/
- https://seandavi.github.io/post/2018/03/extracting-clinical-information-using-the-genomicdatacommons-package/
- The Cancer Data Ecosystem: Data and cloud resources for cancer data science
- The GenomicDataCommons and GEOquery Bioconductor Packages
Note:
- The TCGA data such as TCGA-LUAD are not part of clinical trials (described here).
- Each patient has 4 categories data and the 'case_id' is common to them:
- demographic: gender, race, year_of_birth, year_of_death
- diagnoses: tumor_stage, age_at_diagnosis, tumor_grade
- exposures: cigarettes_per_day, alcohol_history, years_smoked, bmi, alcohol_intensity, weight, height
- main: disease_type, primary_site
- The original download (clinical.tsv file) data contains a column 'treatment_or_therapy' but it has missing values for all patients.
Visualization
GenVisR
ComplexHeatmap
limma
- Differential expression analyses for RNA-sequencing and microarray studies
- Case Study using a Bioconductor R pipeline to analyze RNA-seq data (this is linked from limma package user guide). Here we illustrate how to use two Bioconductor packages - Rsubread' and limma - to perform a complete RNA-seq analysis, including Subread'Bold text read mapping, featureCounts read summarization, voom normalization and limma differential expresssion analysis.
- Unbalanced data, non-normal data, Bartlett's test for equal variance across groups and SAM tests (assumes equal variances just like limma). See this post.
easyRNASeq
Calculates the coverage of high-throughput short-reads against a genome of reference and summarizes it per feature of interest (e.g. exon, gene, transcript). The data can be normalized as 'RPKM' or by the 'DESeq' or 'edgeR' package.
ShortRead
Base classes, functions, and methods for representation of high-throughput, short-read sequencing data.
Rsamtools
The Rsamtools package provides an interface to BAM files.
The main purpose of the Rsamtools package is to import BAM files into R. Rsamtools also provides some facility for file access such as record counting, index file creation, and filtering to create new files containing subsets of the original. An important use case for Rsamtools is as a starting point for creating R objects suitable for a diversity of work flows, e.g., AlignedRead objects in the ShortRead package (for quality assessment and read manipulation), or GAlignments objects in GenomicAlignments package (for RNA-seq and other applications). Those desiring more functionality are encouraged to explore samtools and related software efforts
This package provides an interface to the 'samtools', 'bcftools', and 'tabix' utilities (see 'LICENCE') for manipulating SAM (Sequence Alignment / Map), FASTA, binary variant call (BCF) and compressed indexed tab-delimited (tabix) files.
IRanges
IRanges is a fundamental package (see how many packages depend on it) to other packages like GenomicRanges, GenomicFeatures and GenomicAlignments. The package defines the IRanges class.
The plotRanges() function given in the 'An Introduction to IRanges' vignette shows how to draw an IRanges object.
If we want to make the same plot using the ggplot2 package, we can follow the example in this post. Note that disjointBins() returns a vector the bin number for each bins counting on the y-axis.
flank
The example is obtained from ?IRanges::flank.
ir3 <- IRanges(c(2,5,1), c(3,7,3)) # IRanges of length 3 # start end width # [1] 2 3 2 # [2] 5 7 3 # [3] 1 3 3 flank(ir3, 2) # start end width # [1] 0 1 2 # [2] 3 4 2 # [3] -1 0 2 # Note: by default flank(ir3, 2) = flank(ir3, 2, start = TRUE, both=FALSE) # For example, [2,3] => [2,X] => (..., 0, 1, 2) => [0, 1] # == == flank(ir3, 2, start=FALSE) # start end width # [1] 4 5 2 # [2] 8 9 2 # [3] 4 5 2 # For example, [2,3] => [X,3] => (..., 3, 4, 5) => [4,5] # == == flank(ir3, 2, start=c(FALSE, TRUE, FALSE)) # start end width # [1] 4 5 2 # [2] 3 4 2 # [3] 4 5 2 # Combine the ideas of the previous 2 cases. flank(ir3, c(2, -2, 2)) # start end width # [1] 0 1 2 # [2] 5 6 2 # [3] -1 0 2 # The original statement is the same as flank(ir3, c(2, -2, 2), start=T, both=F) # For example, [5, 7] => [5, X] => ( 5, 6) => [5, 6] # == == flank(ir3, -2, start=F) # start end width # [1] 2 3 2 # [2] 6 7 2 # [3] 2 3 2 # For example, [5, 7] => [X, 7] => (..., 6, 7) => [6, 7] # == == flank(ir3, 2, both = TRUE) # start end width # [1] 0 3 4 # [2] 3 6 4 # [3] -1 2 4 # The original statement is equivalent to flank(ir3, 2, start=T, both=T) # (From the manual) If both = TRUE, extends the flanking region width positions into the range. # The resulting range thus straddles the end point, with width positions on either side. # For example, [2, 3] => [2, X] => (..., 0, 1, 2, 3) => [0, 3] # == # == == == == flank(ir3, 2, start=FALSE, both=TRUE) # start end width # [1] 2 5 4 # [2] 6 9 4 # [3] 2 5 4 # For example, [2, 3] => [X, 3] => (..., 2, 3, 4, 5) => [4, 5] # == # == == == ==
Both IRanges and GenomicRanges packages provide the flank function.
Flanking region is also a common term in High-throughput sequencing. The IGV user guide also has some option related to flanking.
- General tab: Feature flanking regions (base pairs). IGV adds the flank before and after a feature locus when you zoom to a feature, or when you view gene/loci lists in multiple panels.
- Alignments tab: Splice junction track options. The minimum amount of nucleotide coverage required on both sides of a junction for a read to be associated with the junction. This affects the coverage of displayed junctions, and the display of junctions covered only by reads with small flanking regions.
Biostrings
- Manipulate Biological Data Using Biostrings Package Exercises (Part 1)
- Manipulate Biological Data Using Biostrings Package Exercises (Part 2) - it covers global, local & overlap alignments.
GenomicRanges
GenomicRanges depends on IRanges package. See the dependency diagram below.
GenomicFeatues ------- GenomicRanges -+- IRanges -- BioGenomics
| +
+-----+ +- GenomeInfoDb
| |
GenomicAlignments +--- Rsamtools --+-----+
+--- Biostrings
The package defines some classes
- GRanges
- GRangesList
- GAlignments
- SummarizedExperiment: it has the following slots - expData, rowData, colData, and assays. Accessors include assays(), assay(), colData(), expData(), mcols(), ... The mcols() method is defined in the S4Vectors package.
(As of Jan 6, 2015) The introduction in GenomicRanges vignette mentions the GAlignments object created from a 'BAM' file discarding some information such as SEQ field, QNAME field, QUAL, MAPQ and any other information that is not needed in its document. This means that multi-reads don't receive any special treatment. Also pair-end reads will be treated as single-end reads and the pairing information will be lost. This might change in the future.
GenomicAlignments
Counting reads with summarizeOverlaps vignette
library(GenomicAlignments)
library(DESeq)
library(edgeR)
fls <- list.files(system.file("extdata", package="GenomicAlignments"),
recursive=TRUE, pattern="*bam$", full=TRUE)
features <- GRanges(
seqnames = c(rep("chr2L", 4), rep("chr2R", 5), rep("chr3L", 2)),
ranges = IRanges(c(1000, 3000, 4000, 7000, 2000, 3000, 3600, 4000,
7500, 5000, 5400), width=c(rep(500, 3), 600, 900, 500, 300, 900,
300, 500, 500)), "-",
group_id=c(rep("A", 4), rep("B", 5), rep("C", 2)))
features
# GRanges object with 11 ranges and 1 metadata column:
# seqnames ranges strand | group_id
# <Rle> <IRanges> <Rle> | <character>
# [1] chr2L [1000, 1499] - | A
# [2] chr2L [3000, 3499] - | A
# [3] chr2L [4000, 4499] - | A
# [4] chr2L [7000, 7599] - | A
# [5] chr2R [2000, 2899] - | B
# ... ... ... ... ... ...
# [7] chr2R [3600, 3899] - | B
# [8] chr2R [4000, 4899] - | B
# [9] chr2R [7500, 7799] - | B
# [10] chr3L [5000, 5499] - | C
# [11] chr3L [5400, 5899] - | C
# -------
# seqinfo: 3 sequences from an unspecified genome; no seqlengths
olap
# class: SummarizedExperiment
# dim: 11 2
# exptData(0):
# assays(1): counts
# rownames: NULL
# rowData metadata column names(1): group_id
# colnames(2): sm_treated1.bam sm_untreated1.bam
# colData names(0):
assays(olap)$counts
# sm_treated1.bam sm_untreated1.bam
# [1,] 0 0
# [2,] 0 0
# [3,] 0 0
# [4,] 0 0
# [5,] 5 1
# [6,] 5 0
# [7,] 2 0
# [8,] 376 104
# [9,] 0 0
# [10,] 0 0
# [11,] 0 0
Pasilla data. Note that the bam files are not clear where to find them. According to the message, we can download SAM files first and then convert them to BAM files by samtools (Not verify yet).
samtools view -h -o outputFile.bam inputFile.sam
A modified R code that works is
###################################################
### code chunk number 11: gff (eval = FALSE)
###################################################
library(rtracklayer)
fl <- paste0("ftp://ftp.ensembl.org/pub/release-62/",
"gtf/drosophila_melanogaster/",
"Drosophila_melanogaster.BDGP5.25.62.gtf.gz")
gffFile <- file.path(tempdir(), basename(fl))
download.file(fl, gffFile)
gff0 <- import(gffFile, asRangedData=FALSE)
###################################################
### code chunk number 12: gff_parse (eval = FALSE)
###################################################
idx <- mcols(gff0)$source == "protein_coding" &
mcols(gff0)$type == "exon" &
seqnames(gff0) == "4"
gff <- gff0[idx]
## adjust seqnames to match Bam files
seqlevels(gff) <- paste("chr", seqlevels(gff), sep="")
chr4genes <- split(gff, mcols(gff)$gene_id)
###################################################
### code chunk number 12: gff_parse (eval = FALSE)
###################################################
library(GenomicAlignments)
# fls <- c("untreated1_chr4.bam", "untreated3_chr4.bam")
fls <- list.files(system.file("extdata", package="pasillaBamSubset"),
recursive=TRUE, pattern="*bam$", full=TRUE)
path <- system.file("extdata", package="pasillaBamSubset")
bamlst <- BamFileList(fls)
genehits <- summarizeOverlaps(chr4genes, bamlst, mode="Union") # SummarizedExperiment object
assays(genehits)$counts
###################################################
### code chunk number 15: pasilla_exoncountset (eval = FALSE)
###################################################
library(DESeq)
expdata = MIAME(
name="pasilla knockdown",
lab="Genetics and Developmental Biology, University of
Connecticut Health Center",
contact="Dr. Brenton Graveley",
title="modENCODE Drosophila pasilla RNA Binding Protein RNAi
knockdown RNA-Seq Studies",
pubMedIds="20921232",
url="http://www.ncbi.nlm.nih.gov/projects/geo/query/acc.cgi?acc=GSE18508",
abstract="RNA-seq of 3 biological replicates of from the Drosophila
melanogaster S2-DRSC cells that have been RNAi depleted of mRNAs
encoding pasilla, a mRNA binding protein and 4 biological replicates
of the the untreated cell line.")
design <- data.frame(
condition=c("untreated", "untreated"),
replicate=c(1,1),
type=rep("single-read", 2), stringsAsFactors=TRUE)
library(DESeq)
geneCDS <- newCountDataSet(
countData=assay(genehits),
conditions=design)
experimentData(geneCDS) <- expdata
sampleNames(geneCDS) = colnames(genehits)
###################################################
### code chunk number 16: pasilla_genes (eval = FALSE)
###################################################
chr4tx <- split(gff, mcols(gff)$transcript_id)
txhits <- summarizeOverlaps(chr4tx, bamlst)
txCDS <- newCountDataSet(assay(txhits), design)
experimentData(txCDS) <- expdata
We can also check out ?summarizeOverlaps to find some fake examples.
tidybulk
- https://bioconductor.org/packages/release/bioc/html/tidybulk.html
- A Tidy Transcriptomics introduction to RNA-Seq analyses
chromoMap
Rsubread
See this post for about C version of the featureCounts program.
Inference
- How well do RNA-Seq differential gene expression tools perform in a complex eukaryote? A case study in A. thaliana Froussios, Bioinformatics 2019
tximport
- Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences by Soneson, et al. F1000Research. Open peer review.
- http://bioconductor.org/packages/release/bioc/html/tximport.html
$ head -5 quant.sf Name Length EffectiveLength TPM NumReads ENST00000456328.2 1657 1410.79 0.083908 2.46885 ENST00000450305.2 632 410.165 0 0 ENST00000488147.1 1351 1035.94 10.4174 225.073 ENST00000619216.1 68 24 0 0 ENST00000473358.1 712 453.766 0 0
- Another real data from PDMR -> PDX. Select Genomic Analysis -> RNASeq and TPM(genes) column. Consider Patient ID=114348, Specimen ID=004-R, Sample ID=ATHY22, CTEP SDC Code=10038045,
$ head -2 114348_004-R_ATHY22_v2.0.1.4.0_RNASeq.RSEM.genes.results
gene_id transcript_id.s. length effective_length expected_count TPM FPKM
5S_rRNA uc021ofx.1,uc021olu.1,uc021ome.1,uc021oov.1,uc021opl.1,uc021osl.1,uc021ovd.1,uc021ovp.1,uc021pdq.1,uc021pdv.1,uc021pdw.1,uc021ped.1,uc021pff.1,uc021pft.1,uc021phz.1,uc021pmi.1,uc021pmw.1,uc021pnb.1,uc021poq.1,uc021pou.1,uc021ptv.1,uc021pwj.1,uc021pwo.1,uc021pxx.1,uc021qdg.1,uc021qeb.1,uc021qfe.1,uc021qmn.1,uc021qox.1,uc021qqc.1,uc021qqh.1,uc021qsf.1,uc021qvm.1,uc021qzz.1,uc021rcl.1,uc021rij.1,uc021rjt.1,uc021rli.1,uc021rlm.1,uc021rln.1,uc021rmi.1,uc021rmm.1,uc021rva.1,uc021rwx.1,uc021rym.1,uc021sjn.1,uc021snj.1,uc021stz.1,uc021tcu.1,uc021tcv.1,uc021tfq.1,uc021tjk.1,uc021tly.1,uc021tnx.1,uc021tov.1,uc021ubz.1,uc021ujp.1,uc021urx.1,uc021vfh.1,uc021vfo.1,uc021vfu.1,uc021vhy.1,uc021vib.1,uc021vih.1,uc021vjl.1,uc021vjp.1,uc021vju.1,uc021vnc.1,uc021vqb.1,uc021vqp.1,uc021vrv.1,uc021vrx.1,uc021vtp.1,uc021vvw.1,uc021vzi.1,uc021wah.1,uc021wds.1,uc021wek.1,uc021wfa.1,uc021whx.1,uc021wjp.1,uc021wog.1,uc021wpw.1,uc021wuf.1,uc021wuz.1,uc021wyq.1,uc021wzq.1,uc021wzu.1,uc021wzx.1,uc021xai.1,uc021xde.1,uc021xdi.1,uc021xeg.1,uc021xeo.1,uc021xev.1,uc021xez.1,uc021xfx.1,uc021xlz.1,uc021xqf.1,uc021xsq.1,uc021xug.1,uc021xwn.1,uc021xxd.1,uc021xxi.1,uc021xys.1,uc021yks.1,uc021ymh.1,uc021zdo.1,uc021zea.1,uc021zez.1,uc021zfd.1,uc021zft.1,uc021zgc.1,uc021zgn.1,uc022abq.1,uc022akm.1,uc022amt.1,uc022aor.1,uc022atj.1,uc022atr.1,uc022att.1,uc022avo.1,uc022axx.1,uc022bex.1,uc022bia.1,uc022bkv.1,uc022blh.1,uc022bll.1,uc022blx.1,uc022bmb.1,uc022bnt.1,uc022bvb.1,uc022byo.1,uc022bzg.1,uc022cdx.1,uc022cer.1,uc022cgh.1,uc031tkn.1 105.04 2.23 2039.99 49629.87 35353.97
$ R
x = read.delim("114348_004-R_ATHY22_v2.0.1.4.0_RNASeq.RSEM.genes.results")
dim(x)
# [1] 28109 7
names(x)
# [1] "gene_id" "transcript_id.s." "length" "effective_length"
# [5] "expected_count" "TPM" "FPKM"
x[1:3, -2]
# gene_id length effective_length expected_count TPM FPKM
# 1 5S_rRNA 105.04 2.23 2039.99 49629.87 35353.97
# 2 5_8S_rRNA 161.00 21.19 0.00 0.00 0.00
# 3 6M1-18 473.00 302.74 0.00 0.00 0.00
y <- read.delim("114348_004-R_ATHY22_v2.0.1.4.0_RNASeq.RSEM.isoforms.results")
dim(y)
# [1] 78375 8
names(y)
# [1] "transcript_id" "gene_id" "length" "effective_length"
# [5] "expected_count" "TPM" "FPKM" "IsoPct"
y[1:3, -1]
# gene_id length effective_length expected_count TPM FPKM IsoPct
# 1 5S_rRNA 110 3.06 0 0 0 0
# 2 5S_rRNA 133 9.08 0 0 0 0
# 3 5S_rRNA 92 0.00 0 0 0 0
DESeq or edgeR
- DESeq2 method
- DESeq2 with a large number of samples -> use DESEq2 to normalize the data and then use do a Wilcoxon rank-sum test on the normalized counts, for each gene separately, or, even better, use a permutation test. See this post. Or consider the limma-voom method instead, which will handle 1000 samples in a few seconds without the need for extra memory.
- edgeR normalization factor post. Normalization factors are computed using the trimmed mean of M-values (TMM) method; see the paper by Robinson & Oshlack 2010 for more details. Briefly, M-values are defined as the library size-adjusted log-ratio of counts between two libraries. The most extreme 30% of M-values are trimmed away, and the mean of the remaining M-values is computed. This trimmed mean represents the log-normalization factor between the two libraries. The idea is to eliminate systematic differences in the counts between libraries, by assuming that most genes are not DE.
- edgeR From reads to genes to pathways: differential expression analysis of RNA-Seq experiments using Rsubread and the edgeR quasi-likelihood pipeline
- Can I feed TCGA normalized count data to EdgeR?
- counts() function and normalized counts.
- Why use Negative binomial distribution in RNA-Seq data?] and the Presentation by Simon Anders.
- XBSeq2 – a fast and accurate quantification of differential expression and differential polyadenylation. By using simulated datasets, they demonstrated that, overall, XBSeq2 performs equally well as XBSeq in terms of several statistical metrics and both perform better than DESeq2 and edgeR.
- DEBrowser: Interactive Differential Expression Analysis and Visualization Tool for Count Data Alper Kucukural 2018
EBSeq
An R package for gene and isoform differential expression analysis of RNA-seq data
http://www.rna-seqblog.com/analysis-of-ebv-transcription-using-high-throughput-rna-sequencing/
prebs
Probe region expression estimation for RNA-seq data for improved microarray comparability
DEXSeq
Inference of differential exon usage in RNA-Seq
rSeqNP
A non-parametric approach for detecting differential expression and splicing from RNA-Seq data
voomDDA: discovery of diagnostic biomarkers and classification of RNA-seq data
http://www.biosoft.hacettepe.edu.tr/voomDDA/
Pathway analysis
GSOAP
GSOAP: a tool for visualization of gene set over-representation analysis
fgsea: Fast Gene Set Enrichment Analysis
GSEABenchmarkeR: Reproducible GSEA Benchmarking
Towards a gold standard for benchmarking gene set enrichment analysis
hypeR
- hypeR: an R package for geneset enrichment workflows
- Efficient, Scalable, and Reproducible Enrichment Workflows from Bioc2020
GSEPD
GSEPD: a Bioconductor package for RNA-seq gene set enrichment and projection display
BAGSE
BAGSE: a Bayesian hierarchical model approach for gene set enrichment analysis 2020
Pipeline
SARTools
SEQprocess
SEQprocess: a modularized and customizable pipeline framework for NGS processing in R package
pasilla and pasillaBamSubset Data
pasilla - Data package with per-exon and per-gene read counts of RNA-seq samples of Pasilla knock-down by Brooks et al., Genome Research 2011.
pasillaBamSubset - Subset of BAM files untreated1.bam (single-end reads) and untreated3.bam (paired-end reads) from "Pasilla" experiment (Pasilla knock-down by Brooks et al., Genome Research 2011).
BitSeq
Transcript expression inference and differential expression analysis for RNA-seq data. The homepage of Antti Honkela.
ReportingTools
The ReportingTools software package enables users to easily display reports of analysis results generated from sources such as microarray and sequencing data.
sequences
More or less an educational package. It has 2 c and c++ source code. It is used in Advanced R programming and package development.
QuasR
Bioinformatics paper
CRAN packages
ssizeRNA
Sample Size Calculation for RNA-Seq Experimental Design
RnaSeqSampleSize
Shiny app
rbamtools
Provides an interface to functions of the 'SAMtools' C-Library by Heng Li
refGenome
The packge contains functionality for import and managing of downloaded genome annotation Data from Ensembl genome browser (European Bioinformatics Institute) and from UCSC genome browser (University of California, Santa Cruz) and annotation routines for genomic positions and splice site positions.
WhopGenome
Provides very fast access to whole genome, population scale variation data from VCF files and sequence data from FASTA-formatted files. It also reads in alignments from FASTA, Phylip, MAF and other file formats. Provides easy-to-use interfaces to genome annotation from UCSC and Bioconductor and gene ontology data from AmiGO and is capable to read, modify and write PLINK .PED-format pedigree files.
TCGA2STAT
Simple TCGA Data Access for Integrated Statistical Analysis in R
TCGA2STAT depends on Bioconductor package CNTools which cannot be installed automatically.
source("https://bioconductor.org/biocLite.R")
biocLite("CNTools")
install.packages("TCGA2STAT")
The getTCGA() function allows to download various kind of data:
- gene expression which includes mRNA-microarray gene expression data (data.type="mRNA_Array") & RNA-Seq gene expression data (data.type="RNASeq")
- miRNA expression which includes miRNA-array data (data.type="miRNA_Array") & miRNA-Seq data (data.type="miRNASeq")
- mutation data (data.type="Mutation")
- methylation expression (data.type="Methylation")
- copy number changes (data.type="CNA_SNP")
TCGAbiolinks
- https://bioconductor.org/packages/3.12/TCGAbiolinks
- data access through GenomicDataCommons
- The help page of GDCquery does not say it clearly about the option of file.type. How to use TCGAbiolinks to download raw RSEM gene counts for specific type of cancer. file.type= "results" or file.type= "normalized_results".
- An example from Public Data Resources in Bioconductor workshop 2020. According to ?GDCquery, for the legacy data arguments project, data.category, platform and/or file.extension should be used.
library(TCGAbiolinks) library(SummarizedExperiment) query <- GDCquery(project = "TCGA-ACC", data.category = "Gene expression", data.type = "Gene expression quantification", platform = "Illumina HiSeq", file.type = "normalized_results", experimental.strategy = "RNA-Seq", legacy = TRUE) gdcdir <- file.path("Waldron_PublicData", "GDCdata") GDCdownload(query, method = "api", files.per.chunk = 10, directory = gdcdir) # 79 files ACCse <- GDCprepare(query, directory = gdcdir) ACCse class(ACCse) dim(assay(ACCse)) # 19947 x 79 assay(ACCse)[1:3, 1:2] # symbol id length(unique(rownames(assay(ACCse)))) # 19672 rowData(ACCse)[1:2, ] # DataFrame with 2 rows and 3 columns # gene_id entrezgene ensembl_gene_id # <character> <integer> <character> # A1BG A1BG 1 ENSG00000121410 # A2M A2M 2 ENSG00000175899 - HTSeq counts data example
query2 <- GDCquery(project = "TCGA-ACC", data.category = "Transcriptome Profiling", data.type = "Gene Expression Quantification", workflow.type="HTSeq - Counts") gdcdir2 <- file.path("Waldron_PublicData", "GDCdata2") GDCdownload(query2, method = "api", files.per.chunk = 10, directory = gdcdir2) # 79 files ACCse2 <- GDCprepare(query2, directory = gdcdir2) ACCse2 dim(assay(ACCse2)) # 56457 x 79 assay(ACCse2)[1:3, 1:2] # ensembl id rowData(ACCse2)[1:2, ] DataFrame with 2 rows and 3 columns ensembl_gene_id external_gene_name original_ensembl_gene_id <character> <character> <character> ENSG00000000003 ENSG00000000003 TSPAN6 ENSG00000000003.13 ENSG00000000005 ENSG00000000005 TNMD ENSG00000000005.5 - Clinical data
acc_clin <- GDCquery_clinic(project = "TCGA-ACC", type = "Clinical") dim(acc_clin) # [1] 92 71
- TCGAanalyze_DEA(). Differentially Expression Analysis (DEA) Using edgeR Package.
dataNorm <- TCGAbiolinks::TCGAanalyze_Normalization(dataBRCA, geneInfo) dataFilt <- TCGAanalyze_Filtering(tabDF = dataBRCA, method = "quantile", qnt.cut = 0.25) samplesNT <- TCGAquery_SampleTypes(colnames(dataFilt), typesample = c("NT")) samplesTP <- TCGAquery_SampleTypes(colnames(dataFilt), typesample = c("TP")) dataDEGs <- TCGAanalyze_DEA(dataFilt[,samplesNT], dataFilt[,samplesTP],"Normal", "Tumor") # 2nd example dataDEGs <- TCGAanalyze_DEA(mat1 = dataFiltLGG, mat2 = dataFiltGBM, Cond1type = "LGG", Cond2type = "GBM", fdr.cut = 0.01, logFC.cut = 1, method = "glmLRT") - Enrichment analysis
ansEA <– TCGAanalyze_EAcomplete(TFname="DEA genes LGG Vs GBM", RegulonList = rownames(dataDEGs)) TCGAvisualize_EAbarplot(tf = rownames(ansEA$ResBP), GOBPTab = ansEA$ResBP, GOCCTab = ansEA$ResCC, GOMFTab = ansEA$ResMF, PathTab = ansEA$ResPat, nRGTab = rownames(dataDEGs), nBar = 20) - mRNA Analysis Pipeline from GDC documentation.
- Papers
curatedTCGAData
- Curated Data From The Cancer Genome Atlas (TCGA) as MultiAssayExperiment Objects
- data access through ExperimentHub
- provides data from the legacy Firehose pipeline
- Multi-omic Integration and Analysis of cBioPortal and TCGA data with MultiAssayExperiment from Bioc2020
- Public data resources and Bioconductor from Bioc2020
caOmicsV
http://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-016-0989-6 Data from TCGA ws used
Visualize multi-dimentional cancer genomics data including of patient information, gene expressions, DNA methylations, DNA copy number variations, and SNP/mutations in matrix layout or network layout.
Map2NCBI
The GetGeneList() function is useful to download Genomic Features (including gene features/symbols) from NCBI (ftp://ftp.ncbi.nih.gov/genomes/MapView/).
> library(Map2NCBI)
> GeneList = GetGeneList("Homo sapiens", build="ANNOTATION_RELEASE.107", savefiles=TRUE, destfile=path.expand("~/"))
# choose [2], [n], and [1] to filter the build and feature information.
# The destination folder will contain seq_gene.txt, seq_gene.md.gz and GeneList.txt files.
> str(GeneList)
'data.frame': 52157 obs. of 15 variables:
$ tax_id : chr "9606" "9606" "9606" "9606" ...
$ chromosome : chr "1" "1" "1" "1" ...
$ chr_start : num 11874 14362 17369 30366 34611 ...
$ chr_stop : num 14409 29370 17436 30503 36081 ...
$ chr_orient : chr "+" "-" "-" "+" ...
$ contig : chr "NT_077402.3" "NT_077402.3" "NT_077402.3" "NT_077402.3" ...
$ ctg_start : num 1874 4362 7369 20366 24611 ...
$ ctg_stop : num 4409 19370 7436 20503 26081 ...
$ ctg_orient : chr "+" "-" "-" "+" ...
$ feature_name : chr "DDX11L1" "WASH7P" "MIR6859-1" "MIR1302-2" ...
$ feature_id : chr "GeneID:100287102" "GeneID:653635" "GeneID:102466751" "GeneID:100302278" ...
$ feature_type : chr "GENE" "GENE" "GENE" "GENE" ...
$ group_label : chr "GRCh38.p2-Primary" "GRCh38.p2-Primary" "GRCh38.p2-Primary" "GRCh38.p2-Primary" ...
$ transcript : chr "Assembly" "Assembly" "Assembly" "Assembly" ...
$ evidence_code: chr "-" "-" "-" "-" ...
> GeneList$feature_name[grep("^NAP", GeneList$feature_name)]
TCseq: Time course sequencing data analysis
http://bioconductor.org/packages/devel/bioc/html/TCseq.html
UCSCXenaTools
https://github.com/ropensci/UCSCXenaTools
GEO
See the internal link at R-GEO.
GREIN: An interactive web platform for re-analyzing GEO RNA-seq data
GEO2RNAseq
GEO2RNAseq: An easy-to-use R pipeline for complete pre-processing of RNA-seq data
Journals
Biometrical Journal
- https://onlinelibrary.wiley.com/journal/15214036
- Author's Guideline
- Biostatistics (topic) journals from Wiley
Biostatistics
Bioinformatics
Genome Analysis section
BMC Bioinformatics
BioRxiv
PLOS
Software
BRB-SeqTools
https://brb.nci.nih.gov/seqtools/
WebMeV
GeneSpring
RNA-Seq
CCBR Exome Pipeliner
https://ccbr.github.io/Pipeliner/
MOFA: Multi-Omics Factor Analysis
Benchmarking
Essential guidelines for computational method benchmarking
Simulation
Simulate RNA-Seq
- http://en.wikipedia.org/wiki/List_of_RNA-Seq_bioinformatics_tools#RNA-Seq_simulators
- https://popmodels.cancercontrol.cancer.gov/gsr/packages/
Maq
Used by TopHat: discovering splice junctions with RNA-Seq
BEERS/Grant G.R. 2011
http://bioinformatics.oxfordjournals.org/content/27/18/2518.long#sec-2. The simulation method is called BEERS and it was used in the STAR software paper.
For the command line options of <reads_simulator.pl> and more details about the config files that are needed/prepared by BEERS, see this gist.
This can generate paired end data but they are in one FASTA file.
$ sudo apt-get install cpanminus $ sudo cpanm Math::Random $ wget http://cbil.upenn.edu/BEERS/beers.tar $ $ tar -xvf beers.tar # two perl files <make_config_files_for_subset_of_gene_ids.pl> and <reads_simulator.pl> $ $ cd ~/Downloads/ $ mkdir beers_output $ mkdir beers_simulator_refseq && cd "$_" $ wget http://itmat.rum.s3.amazonaws.com/simulator_config_refseq.tar.gz $ tar xzvf simulator_config_refseq.tar.gz $ ls -lth total 1.4G -rw-r--r-- 1 brb brb 44M Sep 16 2010 simulator_config_featurequantifications_refseq -rw-r--r-- 1 brb brb 7.7M Sep 15 2010 simulator_config_geneinfo_refseq -rw-r--r-- 1 brb brb 106M Sep 15 2010 simulator_config_geneseq_refseq -rw-r--r-- 1 brb brb 1.3G Sep 15 2010 simulator_config_intronseq_refseq $ cd ~/Downloads/ $ perl reads_simulator.pl 100 testbeers \ -configstem refseq \ -customcfgdir ~/Downloads/beers_simulator_refseq \ -outdir ~/Downloads/beers_output $ ls -lh beers_output total 3.9M -rw-r--r-- 1 brb brb 1.8K Mar 16 15:25 simulated_reads2genes_testbeers.txt -rw-r--r-- 1 brb brb 1.2M Mar 16 15:25 simulated_reads_indels_testbeers.txt -rw-r--r-- 1 brb brb 1.6K Mar 16 15:25 simulated_reads_junctions-crossed_testbeers.txt -rw-r--r-- 1 brb brb 2.7M Mar 16 15:25 simulated_reads_substitutions_testbeers.txt -rw-r--r-- 1 brb brb 6.3K Mar 16 15:25 simulated_reads_testbeers.bed -rw-r--r-- 1 brb brb 31K Mar 16 15:25 simulated_reads_testbeers.cig -rw-r--r-- 1 brb brb 22K Mar 16 15:25 simulated_reads_testbeers.fa -rw-r--r-- 1 brb brb 584 Mar 16 15:25 simulated_reads_testbeers.log $ wc -l simulated_reads2genes_testbeers.txt 102 simulated_reads2genes_testbeers.txt $ head -4 simulated_reads2genes_testbeers.txt seq.1 GENE.5600 seq.2 GENE.35506 seq.3 GENE.506 seq.4 GENE.34922 $ tail -4 simulated_reads2genes_testbeers.txt seq.97 GENE.4197 seq.98 GENE.8763 seq.99 GENE.19573 seq.100 GENE.18830 $ wc -l simulated_reads_indels_testbeers.txt 36131 simulated_reads_indels_testbeers.txt $ head -2 simulated_reads_indels_testbeers.txt chr1:6052304-6052531 25 1 G chr2:73899436-73899622 141 3 ATA $ tail -2 simulated_reads_indels_testbeers.txt chr4:68619532-68621804 1298 -2 AA chr21:32554738-32554962 174 1 T $ wc -l simulated_reads_substitutions_testbeers.txt 71678 simulated_reads_substitutions_testbeers.txt $ head -2 simulated_reads_substitutions_testbeers.txt chr22:50902963-50903167 50903077 G->A chr1:6052304-6052531 6052330 G->C $ wc -l simulated_reads_junctions-crossed_testbeers.txt 49 simulated_reads_junctions-crossed_testbeers.txt $ head -2 simulated_reads_junctions-crossed_testbeers.txt seq.1a chrX:49084601-49084713 seq.1b chrX:49084909-49086682 $ cat beers_output/simulated_reads_testbeers.log Simulator run: 'testbeers' started: Thu Mar 16 15:25:39 EDT 2017 num reads: 100 readlength: 100 substitution frequency: 0.001 indel frequency: 0.0005 base error: 0.005 low quality tail length: 10 percent of tails that are low quality: 0 quality of low qulaity tails: 0.8 percent of alt splice forms: 0.2 number of alt splice forms per gene: 2 stem: refseq sum of gene counts: 3,886,863,063 sum of intron counts = 1,304,815,198 sum of intron counts = 2,365,472,596 intron frequency: 0.355507598105262 padded intron frequency: 0.52453796437909 finished at Thu Mar 16 15:25:58 EDT 2017 $ wc -l simulated_reads_testbeers.fa 400 simulated_reads_testbeers.fa $ head simulated_reads_testbeers.fa >seq.1a CGAAGAAGGACCCAAAGATGACAAGGCTCACAAAGTACACCCAGGGCAGTTCATACCCCATGGCATCTTGCATCCAGTAGAGCACATCGGTCCAGCCTTC >seq.1b GCTCGAGCTGTTCCTTGGACGAATGCACAAGACGTGCTACTTCCTGGGATCCGACATGGAAGCGGAGGAGGACCCATCGCCCTGTGCATCTTCGGGATCA >seq.2a GCCCCAGCAGAGCCGGGTAAAGATCAGGAGGGTTAGAAAAAATCAGCGCTTCCTCTTCCTCCAAGGCAGCCAGACTCTTTAACAGGTCCGGAGGAAGCAG >seq.2b ATGAAGCCTTTTCCCATGGAGCCATATAACCATAATCCCTCAGAAGTCAAGGTCCCAGAATTCTACTGGGATTCTTCCTACAGCATGGCTGATAACAGAT >seq.3a CCCCAGAGGAGCGCCACCTGTCCAAGATGCAGCAGAACGGCTACGAAAATCCAACCTACAAGTTCTTTGAGCAGATGCAGAACTAGACCCCCGCCACAGC # Take a look at the true coordinates $ head -4 simulated_reads_testbeers.bed # one-based coords and contains both endpoints of each span chrX 49084529 49084601 + chrX 49084713 49084739 + chrX 49084863 49084909 + chrX 49086682 49086734 + $ head -4 simulated_reads_testbeers.cig # has a cigar string representation of the mapping coordinates, and a more human readable representation of the coordinates seq.1a chrX 49084529 73M111N27M 49084529-49084601, 49084713-49084739 + CGAAGAAGGACCCAAAGATGACAAGGCTCACAAAGTACACCCAGGGCAGTTCATACCCCATGGCATCTTGCATCCAGTAGAGCACATCGGTCCAGCCTTC seq.1b chrX 49084863 47M1772N53M 49084863-49084909, 49086682-49086734 - GCTCGAGCTGTTCCTTGGACGAATGCACAAGACGTGCTACTTCCTGGGATCCGACATGGAAGCGGAGGAGGACCCATCGCCCTGTGCATCTTCGGGATCA seq.2a chr1 183516256 100M 183516256-183516355 - GCCCCAGCAGAGCCGGGTAAAGATCAGGAGGGTTAGAAAAAATCAGCGCTTCCTCTTCCTCCAAGGCAGCCAGACTCTTTAACAGGTCCGGAGGAAGCAG seq.2b chr1 183515275 100M 183515275-183515374 + ATGAAGCCTTTTCCCATGGAGCCATATAACCATAATCCCTCAGAAGTCAAGGTCCCAGAATTCTACTGGGATTCTTCCTACAGCATGGCTGATAACAGAT $ wc -l simulated_reads_testbeers.fa 400 simulated_reads_testbeers.fa $ wc -l simulated_reads_testbeers.bed 247 simulated_reads_testbeers.bed $ wc -l simulated_reads_testbeers.cig 200 simulated_reads_testbeers.cig
Flux Sammeth 2010
SimNGS
SimSeq
A data-based simulation algorithm for rna-seq data. The vector of read counts simulated for a given experimental unit has a joint distribution that closely matches the distribution of a source rna-seq dataset provided by the user.
empiricalFDR.DESeq2
http://biorxiv.org/content/early/2014/12/05/012211
The key function is simulateCounts, which takes a fitted DESeq2 data object as an input and returns a simulated data object (DESeq2 class) with the same sample size factors, total counts and dispersions for each gene as in real data, but without the effect of predictor variables.
Functions fdrTable, fdrBiCurve and empiricalFDR compare the DESeq2 results obtained for the real and simulated data, compute the empirical false discovery rate (the ratio of the number of differentially expressed genes detected in the simulated data and their number in the real data) and plot the results.
polyester
http://biorxiv.org/content/early/2014/12/05/012211
Given a set of annotated transcripts, polyester will simulate the steps of an RNA-seq experiment (fragmentation, reverse-complementing, and sequencing) and produce files containing simulated RNA-seq reads.
Input: reference FASTA file (containing names and sequences of transcripts from which reads should be simulated) OR a GTF file denoting transcript structures, along with one FASTA file of the DNA sequence for each chromosome in the GTF file.
Output: FASTA files. Reads in the FASTA file will be labeled with the transcript from which they were simulated.
Too many dependencies. Got an error in installation.. It seems it has not considered splice junctions.
seqgendiff
Data-based RNA-seq simulations by binomial thinning
RSEM
- Cbioportal -> Breast Invasive Carcinoma (TCGA, PanCancer Atlas) -> Explore selected study -> Original data -> PanCanAtlas Publications -> RNA (Final) - EBPlusPlusAdjustPANCAN_IlluminaHiSeq_RNASeqV2.geneExp.tsv
Simulate DNA-Seq
- Software list - https://popmodels.cancercontrol.cancer.gov/gsr/packages/
- A comparison of tools for the simulation of genomic next-generation sequencing data Merly Escalona 2016
wgsim
- Used by Cleaning clinical genomic data: Simple identification and removal of recurrently miscalled variants in single genomes bioRxiv 2017
- (How to) Simulate reads using a reference genome ALT contig
- [http://research.cs.wisc.edu/wham/comparison-using-wgsim/ Comparing WHAM with BWA using wgsim
- http://biobits.org/samtools_primer.html
dwgssim
- Whole Genome Simulator for Next-Generation Sequencing
- Example: How to generate variant calls with bcftools
NEAT
- https://github.com/zstephens/neat-genreads
- Simulating next-generation sequencing datasets from empirical mutation and sequencing models Zachary Stephens, 2016
- If I set 10 as the coverage rate and read length 101, the generated fq file is about 34GB (3.3GB * 10) for each one of the pairs.
DNA aligner accuracy: BWA, Bowtie, Soap and SubRead tested with simulated reads
http://genomespot.blogspot.com/2014/11/dna-aligner-accuracy-bwa-bowtie-soap.html
$ head simDNA_100bp_16del.fasta >Pt-0-100 TGGCGAACGCGGGAATTGACCGCGATGGTGATTCACATCACTCCTAATCCACTTGCTAATCGCCCTACGCTACTATCATTCTTT >Pt-10-110 GCGGGATTGAACCCGATTGAATTCCAATCACTGCTTAATCCACTTGCTACATCGCCCTACGTACTATCTATTTTTTTGTATTTC >Pt-20-120 GAACCCGCGATGAATTCAATCCACTGCTACCATTGGCTACATCCGCCCCTACGCTACTCTTCTTTTTTGTATGTCTAAAAAAAA >Pt-30-130 TGGTGAATCACAATCACTGCCTAACCATTGGCTACATCCGCCCCTACGCTACACTATTTTTTGTATTGCTAAAAAAAAAAATAA >Pt-40-140 ACAACACTGCCTAATCCACTTGGCTACTCCGCCCCTAGCTACTATCTTTTTTTGTATTTCTAAAAAAAAAAAATCAATTTCAAT
Simulate Whole genome
- DWGSIM mentioned by Variant Callers for Next-Generation Sequencing Data: A Comparison Study. For its usage, see http://davetang.org/wiki/tiki-index.php?page=DWGSIM
Simulate whole exome
- https://www.biostars.org/p/66714/ (no final answer)
- Wessim: a whole-exome sequencing simulator based on in silico exome capture Sangwoo Kim 2013 & software
SCSIM
SCSIM: Jointly simulating correlated single-cell and bulk next-generation DNA sequencing data
Variant simulator
Convert FASTA to FASTQ
It is interesting to note that the simulated/generated FASTA files can be used by alignment/mapping tools like BWA just like FASTQ files.
If we want to convert FASTA files to FASTQ files, use https://code.google.com/archive/p/fasta-to-fastq/. The quality score 'I' means 40 (the highest) by Sanger (range [0,40]). See https://en.wikipedia.org/wiki/FASTQ_format. The Wikipedia website also mentions FASTQ read simulation tools and a comparison of these tools.
$ cat test.fasta >Pt-0-50 TGGCGAACGACGGGAATACCCGGAGGTGAATTCAAATCCACT >Pt-10-60 GACGGAATTGAACCCGATGGGATACAATCCACTGCCTTATCC >Pt-20-70 GAACCCGCGATGGTGTCACAATCCACTCTTAACCATTGCTAC >Pt-30-80 GGTGAATTCACAATCCACTGCCTTACCACTTGGCTACCCCCT >Pt-40-90 AATCCACTGCCTTATCCACTGGCTACATCCCTACGCTACTAT $ perl ~/Downloads/fasta_to_fastq.pl test.fasta @Pt-0-50 TGGCGAACGACGGGAATACCCGGAGGTGAATTCAAATCCACT + IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII @Pt-10-60 GACGGAATTGAACCCGATGGGATACAATCCACTGCCTTATCC + IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII @Pt-20-70 GAACCCGCGATGGTGTCACAATCCACTCTTAACCATTGCTAC + IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII @Pt-30-80 GGTGAATTCACAATCCACTGCCTTACCACTTGGCTACCCCCT + IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII @Pt-40-90 AATCCACTGCCTTATCCACTGGCTACATCCCTACGCTACTAT + IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
Alternatively we can use just one line of code by awk
$ awk 'BEGIN {RS = ">" ; FS = "\n"} NR > 1 {print "@"$1"\n"$2"\n+"; for(c=0;c<length($2);c++) printf "H"; printf "\n"}' \
test.fasta > test.fq
$ cat test.fq
@Pt-0-50
TGGCGAACGACGGGAATACCCGGAGGTGAATTCAAATCCACT
+
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
@Pt-10-60
GACGGAATTGAACCCGATGGGATACAATCCACTGCCTTATCC
+
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
@Pt-20-70
GAACCCGCGATGGTGTCACAATCCACTCTTAACCATTGCTAC
+
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
@Pt-30-80
GGTGAATTCACAATCCACTGCCTTACCACTTGGCTACCCCCT
+
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
@Pt-40-90
AATCCACTGCCTTATCCACTGGCTACATCCCTACGCTACTAT
+
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
Change the 'H' to the quality score value that you need (Depending what phred score scale you are using).
Simulate genetic data
‘Simulating genetic data with R: an example with deleterious variants (and a pun)’
PDX/Xenograft
- https://en.wikipedia.org/wiki/Patient_derived_xenograft
- Are special read alignment strategies necessary and cost-effective when handling sequencing reads from patient-derived tumor xenografts? by Tso et al, BMC Genomics, 2014.
- Map xenograft reads by Array Suite
- PDXFinder includes PDMR as one of 6 providers. The website source code https://github.com/PDXFinder/pdxfinder.
- A Colorectal Carcinoma example contains 'genomic data'. Each row represents one seq position. No raw FASTQ files available.
- NCI Patient-Derived Models Repository (PDMR).
- ftp://dctdftp.nci.nih.gov/pub/pdm/
- The ftp link can be obtained by clicking 'PDMR Models' -> 'PDMR database' -> Click here to access the PDMR Database -> 'Genomic Analysis'.
- There will be three tabs: NCI Cancer Genome Panel (4246 records), Whole Exome Sequence (785 records) and RNASeq (807 records).
- For Whole Exome Sequence, VCF was provided. For RNASeq, TPM files per genes or isoforms are available.
- The RNASeq Transcriptome Data Analysis Pipeline and Specifications and Whole Exome Sequencing Data Analysis Pipeline and Specifications are available under SOPs. RNASeq TPM data.
- Reproducible Bioinformatics Project
- Xenome a tool for classifying reads from xenograft samples, Thomas Conway et al 2012.
#!/bin/bash module load gossamer xenome index -M 24 -T 16 -P idx \ -H $HOME/igenomes/Mus_musculus/UCSC/mm9/Sequence/WholeGenomeFasta/genome.fa \ -G $HOME/igenomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa
- Next-Generation Sequence Analysis of Cancer Xenograft Models by Fernando J. Rossello et al 2013.
- Whole transcriptome profiling of patient-derived xenograft models as a tool to identify both tumor and stromal specific biomarkers Bradford et al 2016
- An open-source application for disambiguating two species in next generation sequencing data from grafted samples by Ahdesmäki MJ et al 2016. Disambiguation
- Next-Generation Sequencing Analysis and Algorithms for PDX and CDX Models by Garima Khandelwal et al 2017. bamcmp software.
- Computational approach to discriminate human and mouse sequences in patient-derived tumour xenografts by Maurizio Callari et al 2018. Both RNA-Seq and DNA-Seq are considered. Software ICRG.
- XenofilteR: computational deconvolution of mouse and human reads in tumor xenograft sequence data Kluin et al 2018. Software in github.
- Whole genomes define concordance of matched primary, xenograft, and organoid models of pancreas cancer Gendoo et al 2019
- Integrative Pharmacogenomics Analysis of Patient-Derived Xenografts 2019
RNA-Seq
- Platform GPL16791 Illumina HiSeq 2500
- GPL25431 Illumina HiSeq 4000 (Homo sapiens; Mus musculus)
- GSE118197
- GPL18573 Illumina NextSeq 500
- GSE106336 MYC regulates ductal-neuroendocrine lineage plasticity in pancreatic ductal adenocarcinoma associated with poor outcome and chemoresistance
- Mastering RNA-Seq (NGS Data Analysis) - A Critical Approach To Transcriptomic Data Analysis. This is posted on rna-seqblog. See the link Whole transcriptome profiling of patient-derived xenograft models as a tool to identify both tumor and stromal specific biomarkers.
- PDXGEM: patient-derived tumor xenograft-based gene expression model for predicting clinical response to anticancer therapy in cancer patients
DNA-Seq
- Endocrine-Therapy-Resistant ESR1 Variants Revealed by Genomic Characterization of Breast-Cancer-Derived Xenografts
- https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=studies&f=study&term=xenograft&go=Go
- ERP021871
MAF (TCGA, GDC)
- See D3Oncoprint manual. D3oncoprint uses java to create the interface (java -jar D3Oncoprint.jar). The result is shown in a browser. The top part is a heatmap and the bottom part is a table. The heatmap is good only for a small number of variants and there is no way to zoom in if we have a lot of variants (eg 1000). The phenotypes file is optional. The following 3 columns are required and the others are for tooltip and table. The column name from a VCF file input is given below.
- Gene column: GENE
- Variant type column: ExonicFunc.refGene
- Aminoacid change column: AAChange.refGene (this defines the name of variant. The variant uses the single letter aminoacid codes. This information is only used in HTML table. For example if AAChange.refGene=UGT1A7:NM_019077:exon1:c.T387G:p.N129K, then the variant will be N129K. My observation is the variant name for variant type 'frameshift_deletion' is altered in table; R128fs -> R128Xfs)
- https://www.reddit.com/r/bioinformatics/comments/cmdza5/vcf_and_maf_file_formats/ which links to
- The Standing Operating Procedure (SOP) describes procedures for generating consensus/intersect variants calls for reporting in the NCI Patient-Derived Models database
- Convert vcf to maf. https://hpc.nih.gov/apps/vcf2maf.html
- maftools : Summarize, Analyze and Visualize MAF Files (oncoplots)
- https://www.bioconductor.org/packages/release/bioc/html/TCGAbiolinks.html
- Use TCGAbiolinks package to download maf file.
DNA Seq Data
NIH
- Go to SRA/Sequence Read Archiveand type the keywords 'Whole Genome Sequencing human'. An example of the procedures to search whole genome sequencing data from human samples:
- Enter 'Whole Genome Sequencing human' in ncbi/sra search sra objects at http://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=search_obj
- The webpage will return the result in terms of SRA experiments, SRA studies, Biosamples, GEO datasets. I pick SRA studies from Public Access.
- The result is sorted by the Accession number (does not take the first 3 letters like DRP into account). The Accession number has a format SRPxxxx. So I just go to the Last page (page 98)
- I pick the first one Accession:SRP066837 from this page. The page shows the Study type is whole genome sequence. http://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?study=SRP066837
- (Important trick) Click the number next to Run. It will show a summary (SRR #, library name, MBases, age, biomaterial provider, isolate and sex) about all samples. http://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP066837
- Download the raw data from any one of them (eg SRR2968056). For whole genome, the Strategy is WGS. For whole exome, the Strategy is called WXS.
- Search the keywords 'nonsynonymous' and 'human' in PMC
Use SRAToolKit instead of wget to download
Don't use the wget command since it requires the specification of right http address.
Downloading SRA data using command line utilities
SRA2R - a package to import SRA data directly into R.
(Method 1) Use the fastq-dump command. For example, the following command (modified from the document will download the first 5 reads and save it to a file called <SRR390728.fastq> (NOT sra format) in the current directory.
/opt/RNA-Seq/bin/sratoolkit.2.3.5-2-ubuntu64/bin/fastq-dump -X 5 SRR390728 -O . # OR /opt/RNA-Seq/bin/sratoolkit.2.3.5-2-ubuntu64/bin/fastq-dump --split-3 SRR390728 # no progress bar
This will download the files in FASTQ format.
(Method 2) If we need to downloading by wget or FTP (works for ‘SRR’, ‘ERR’, or ‘DRR’ series):
wget ftp://ftp-trace.ncbi.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR304/SRR304976/SRR304976.sra
It will download the file in SRA format. In the case of SRR590795, the sra is 240M and fastq files are 615*2MB.
(Method 3) Download Ubuntu x86_64 tarball from http://downloads.asperasoft.com/en/downloads/8?list
brb@T3600 ~/Downloads $ tar xzvf aspera-connect-3.6.2.117442-linux-64.tar.gz aspera-connect-3.6.2.117442-linux-64.sh brb@T3600 ~/Downloads $ ./aspera-connect-3.6.2.117442-linux-64.sh Installing Aspera Connect Deploying Aspera Connect (/home/brb/.aspera/connect) for the current user only. Restart firefox manually to load the Aspera Connect plug-in Install complete. brb@T3600 ~/Downloads $ ~/.aspera/connect/bin/ascp -QT -l640M \ -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh \ [email protected]:/sra/sra-instant/reads/ByRun/sra/SRR/SRR590/SRR590795/SRR590795.sra . SRR590795.sra 100% 239MB 535Mb/s 00:06 Completed: 245535K bytes transferred in 7 seconds (272848K bits/sec), in 1 file. brb@T3600 ~/Downloads $
Aspera is typically 10 times faster than FTP according to the website. For this case, wget takes 12s while ascp uses 7s.
Note that the URL on the website's is wrong. I got the correct URL from emailing to ncbi help. Google: ascp "[email protected]"
SRAdb package
https://bioconductor.org/packages/release/bioc/html/SRAdb.html
First we install some required package for XML and RCurl.
sudo apt-get update sudo apt-get install libxml2-dev sudo apt-get install libcurl4-openssl-dev
and then
source("https://bioconductor.org/biocLite.R")
biocLite("SRAdb")
SRA
Only the cancer types with expected cases > 10^5 in the US in 2015 are considered here. http://www.cancer.gov/types/common-cancers
SRA Explorer
SRP056969
- Inference of RNA decay rate from transcriptional profiling highlights the regulatory programs of Alzheimer’s disease
- RNA-Seq reveals mRNA stability a marker in Alzheimer’s patients
- REMBRANDTS: REMoving Bias from Rna-seq ANalysis of Differential Transcript Stability
SRP066363 - lung cancer
- Platform: GPL11154 Illumina HiSeq 2000 (Homo sapiens)
- Overall design: RNAseq and DNA copy number analysis of H1975 cells
- Strategy: 6 RNA-Seq and 3 Whole exome. Paired. 9 samples
- http://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP066363
- http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE74866
SRP015769 or SRP062882 - prostate cancer
- http://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP015769 5 are from normal and 5 are from tumor. Whole Exome Seq.
- http://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP062882 6 normal and the rest are tumor.
SRP053134 - breast cancer
Look at the MBases value column. It determines the coverage for each run.
SRP050992 single cell RNA-Seq
Used in Design and computational analysis of single-cell RNA-sequencing experiments
Single cell RNA-Seq
- Exploiting single-cell expression to characterize co-expression replicability
- NGS -> Single cell RNA-Seq
SRP040626 or SRP040540 - Colon and rectal cancer
- http://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP040626
- http://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP040540
OmicIDX
Tutorials
See the BWA section.
Whole Exome Seq
- 1000genomes. 1000genomes and tcga are two places to get vcf files too.
- Review of Current Methods, Applications, and Data Management for the Bioinformatics Analysis of Whole Exome Sequencing (Bao 2014)
- A framework for variation discovery and genotyping using next-generation DNA sequencing data. See the table 1 there.
- Some data from SRA repository.
Whole Genome Seq
- BioProject and
- Search: filter by 'Homo sapiens wgs'
- Project data type: Genome sequencing
- Click 'Date' to sort by it
- 20 hits as of 1/11/2017 (many of them do not have data)
- PRJNA352450 18 experiments
- SRX2341045 20.5M spots, 4.1G bases, 1.8Gb
- PRJNA343545 48 experiments
- SRX2187657 417.4M spots, 126.1G bases, 55.1Gb
- 309109 5 experiments
- SRX1538498 1.7M spots, 346.1M bases, 99.7Mb
- PRJNA289286 5 experiments
- SRX1100298 504M spots, 101.8G bases, 45.3Gb
- PRJNA260389 27 experiments
- SRX699196 34,099,675 spots, 6.8G bases, 4.3Gb size
- 248553 3 experiments
- SRX1026041 1.2G spots, 250.9G bases, 114.2Gb
- 210123 26 experiments
- SRX318496 173.7M spots, 34.7G bases, 23Gb
- 43433 3 experiments (ABI SOLiD System 3.0)
- SRX017230 851.1M spots, 85.1G bases, 68.8Gb
SraRunTable.txt
- http://www.ncbi.nlm.nih.gov/sra/?term=SRA059511
- http://www.ncbi.nlm.nih.gov/sra/SRX194938[accn] and click SRP004077
- http://trace.ncbi.nlm.nih.gov/Traces/sra/?study=SRP004077 and click Runs from the RHS
- http://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP004077 and click RunInfoTable
Note that (For this study, it has 2377 rows)
- Column A (AssemblyName_s) eg GRCh37
- Column I (library_name_s) eg
- column N (header=Run_s) shows all SRR or ERR accession numbers.
- Column P (Sample_Name)
- Column Y (header=Assay_Type_s) shows WGS.
- Column AB (LibraryLayout_s): PAIRED
Public Data
Public data resources and Bioconductor from Bioc2020.
| Package name | Object class |
|---|---|
| GEOquery | SummarizedExperiment |
| GenomicDataCommons | GDCQuery |
| TCGAbiolinks curatedTCGAData |
RangedSummarizedExperiment MultiAssayExperimentObjects |
| recount | RangedSummarizedExperiment |
| curatedMetagenomicData | ExperimentHub |
ISB Cancer Genomics Cloud (ISB-CGC)
https://isb-cgc.appspot.com/ Leveraging Google Cloud Platform for TCGA Analysis
The ISB Cancer Genomics Cloud (ISB-CGC) is democratizing access to NCI Cancer Data (TCGA, TARGET, CCLE) and coupling it with unprecedented computational power to allow researchers to explore and analyze this vast data-space.
CCLE
Next-generation characterization of the Cancer Cell Line Encyclopedia 2019
It has 1000+ cell lines profiled with different -omics including DNA methylation, RNA splicing, as well as some proteomics (and lots more!).ß
NCI's Genomic Data Commons (GDC)/TCGA
The GDC supports several cancer genome programs at the NCI Center for Cancer Genomics (CCG), including The Cancer Genome Atlas (TCGA), Therapeutically Applicable Research to Generate Effective Treatments (TARGET), and the Cancer Genome Characterization Initiative (CGCI).
- NCI's GDC - Genomic Data Commons Data Portal. Researchers can access over 3 PB of bigData from projects like CPTAC, TARGET and of course TCGA.
- MDAnderson download, FAQ, Available Software
- Working with TCGAbiolinks package
- edgeR is used to find DE; see the TCGAanalyze_DEA() function and TCGAanalyze: Analyze data from TCGA
- GEO accession data GSE62944 as a SummarizedExperiment and GEO website.
- BioXpress: an integrated RNA-seq-derived gene expression database for pan-cancer analysis.
- https://github.com/srp33/TCGA_RNASeq_Clinical
- MOGSA: integrative single sample gene-set analysis of multiple omics data 2019. The data was obtained from TCGA and NCI60.
- RNA Sequencing of the NCI-60: Integration into CellMiner and CellMiner CDB
- Integrating Multidimensional Data for Clustering Analysis With Applications to Cancer Patient Data 2020
GenomicDataCommons package
- See here
- GenomicDataCommons Example: UUID to TCGA and TARGET Barcode Translation
Case studies
Expression of the SARS-CoV-2 cell receptor gene ACE2 in a wide variety of human tissues
GTEx
- Genotype-Tissue Expression (GTEx) project
- recount workflow: accessing over 70,000 human RNA-seq samples with Bioconductor
- Basal Contamination of Bulk Sequencing: Lessons from the GTEx dataset
DepMap
- Dependency Map (DepMap) portal is to empower the research community to make discoveries related to cancer vulnerabilities by providing open access to key cancer dependencies analytical and visualization tools.
- depmap: Cancer Dependency Map Data Package
Sharing data
- NCI Data Sharing
- Ten quick tips for sharing open genomic data Brown et al, PLOS 2018
Gene set analysis
Hypergeometric test
- A course from bioinformatics.ca and over-represented pathway.
- How informative are enrichment analyses really?
Next-generation sequencing data
Forums
- Biostars source code https://github.com/ialbert/biostar-central
- https://support.bioconductor.org/ (powered by Biostar too)
Misc
Advice
- Ten great papers for biologists starting out in computational biology
- Bioinformatics advice I wish I learned 10 years ago from NIH
High Performance
Cloud Computing
- Falco: A quick and flexible single-cell RNA-seq processing framework on the cloud
- Getting started with Bioconductor in the cloud
- miCloud: a bioinformatics cloud for seamless execution of complex NGS data analysis pipelines
Merge different datasets (different genechips)
Normalization
- How data analysis affects power, reproducibility and biological insight of RNA-seq studies in complex datasets 2015
- Comparing the normalization methods for the differential analysis of Illumina high-throughput RNA-Seq data 2015
- Expression analysis of RNA sequencing data from human neural and glial cell lines depends on technical replication and normalization methods 2018
- A statistical normalization method and differential expression analysis for RNA-seq data between different species Zhou 2019
- RNA-seq analysis some notes from Ming Tang
- A protocol to evaluate RNA sequencing normalization methods Abrams et al 2019. It concluded that transcripts per million (TPM) was the best performing normalization method based on its preservation of biological signal as compared to the other methods tested.
- Data pre-processing 6.6.3 Normalising gene expression distributions. log2 CPM (counts per million) was considered. edgeR::calcNormFactors(x, method = "TMM") was used to normalize counts.
- CrossNorm: a novel normalization strategy for microarray data in cancers Cheng 2016
Ensembl
- http://useast.ensembl.org/index.html
- Training
- Ensembl to gene symbol using ensembl_to_symbol().
- A Tidy Transcriptomics introduction to RNA sequencing analyses (Bioc 2020).
- "org.Hs.eg.db" package
- "biomaRt" package
Ensembl is a genome browser for vertebrate genomes that supports research in comparative genomics, evolution, sequence variation and transcriptional regulation. Ensembl annotate genes, computes multiple alignments, predicts regulatory function and collects disease data. Ensembl tools include BLAST, BLAT, BioMart and the Variant Effect Predictor (VEP) for all supported species.
How to use UCSC Table Browser
- An instruction from BitSeq software
- How To Get Bed File Containing Exons Of Canonical Transcripts And Their Corresponding Gene Symbols
- Where to download refseq gene coding regions data?
- http://genome.ucsc.edu/cgi-bin/hgTables OR
- Download refGene.txt.gz file from UCSC directly using http links
- Where To Download Genome Annotation Including Exon, Intron, Utr, Intergenic Information?
File:Tablebrowser.png File:Tablebrowser2.png
Note
- the UCSC browser will return the output on browser by default. Users need to use the browser to save the file with self-chosen file name.
- the output does not have a header
- The bed format is explained in https://genome.ucsc.edu/FAQ/FAQformat.html#format1
If I select "Whole Genome", I will get a file with 75,893 rows. If I choose "Coding Exons", I will get a file with 577,387 rows.
$ wc -l hg38Tables.bed 75893 hg38Tables.bed $ head -2 hg38Tables.bed chr1 67092175 67134971 NM_001276352 0 - 67093579 67127240 0 9 1429,70,145,68,113,158,92,86,42, 0,4076,11062,19401,23176,33576,34990,38966,42754, chr1 201283451 201332993 NM_000299 0 + 201283702 201328836 0 15 453,104,395,145,208,178,63,115,156,177,154,187,85,107,2920, 0,10490,29714,33101,34120,35166,36364,36815,38526,39561,40976,41489,42302,45310,46622, $ tail -2 hg38Tables.bed chr22_KI270734v1_random 131493 137393 NM_005675 0 + 131645 136994 0 5 262,161,101,141,549, 0,342,3949,4665,5351, chr22_KI270734v1_random 138078 161852 NM_016335 0 - 138479 161586 0 15 589,89,99,176,147,93,82,80,117,65,150,35,209,313,164, 0,664,4115,5535,6670,6925,8561,9545,10037,10335,12271,12908,18210,23235,23610, $ wc -l hg38CodingExon.bed 577387 hg38CodingExon.bed $ head -2 hg38CodingExon.bed chr1 67093579 67093604 NM_001276352_cds_0_0_chr1_67093580_r 0 - chr1 67096251 67096321 NM_001276352_cds_1_0_chr1_67096252_r 0 - $ tail -2 hg38CodingExon.bed chr22_KI270734v1_random 156288 156497 NM_016335_cds_12_0_chr22_KI270734v1_random_156289_r 0 - chr22_KI270734v1_random 161313 161586 NM_016335_cds_13_0_chr22_KI270734v1_random_161314_r 0 - # Focus on one NCBI refseq (https://www.ncbi.nlm.nih.gov/nuccore/444741698) $ grep NM_001276352 hg38Tables.bed chr1 67092175 67134971 NM_001276352 0 - 67093579 67127240 0 9 1429,70,145,68,113,158,92,86,42, 0,4076,11062,19401,23176,33576,34990,38966,42754, $ grep NM_001276352 hg38CodingExon.bed chr1 67093579 67093604 NM_001276352_cds_0_0_chr1_67093580_r 0 - chr1 67096251 67096321 NM_001276352_cds_1_0_chr1_67096252_r 0 - chr1 67103237 67103382 NM_001276352_cds_2_0_chr1_67103238_r 0 - chr1 67111576 67111644 NM_001276352_cds_3_0_chr1_67111577_r 0 - chr1 67115351 67115464 NM_001276352_cds_4_0_chr1_67115352_r 0 - chr1 67125751 67125909 NM_001276352_cds_5_0_chr1_67125752_r 0 - chr1 67127165 67127240 NM_001276352_cds_6_0_chr1_67127166_r 0 -
This can be compared to refGene(?) directly downloaded via http
$ wget -c -O hg38.refGene.txt.gz http://hgdownload.soe.ucsc.edu/goldenPath/hg38/database/refGene.txt.gz --2018-10-09 15:44:43-- http://hgdownload.soe.ucsc.edu/goldenPath/hg38/database/refGene.txt.gz Resolving hgdownload.soe.ucsc.edu (hgdownload.soe.ucsc.edu)... 128.114.119.163 Connecting to hgdownload.soe.ucsc.edu (hgdownload.soe.ucsc.edu)|128.114.119.163|:80... connected. HTTP request sent, awaiting response... 200 OK Length: 7457957 (7.1M) [application/x-gzip] Saving to: ‘hg38.refGene.txt.gz’ hg38.refGene.txt.gz 100%[===============================================================>] 7.11M 901KB/s in 10s 2018-10-09 15:44:54 (708 KB/s) - ‘hg38.refGene.txt.gz’ saved [7457957/7457957] $ zcat hg38.refGene.txt.gz | wc -l 75893 15:45PM /tmp$ zcat hg38.refGene.txt.gz | head -2 1072 NM_003288 chr20 + 63865227 63891545 63865365 63889945 7 63865227,63869295,63873667,63875815,63882718,63889189,63889849, 63865384,63869441,63873816,63875875,63882820,63889238,63891545, 0 TPD52L2 cmpl cmpl 0,1,0,2,2,2,0, 1815 NR_110164 chr2 + 161244738 161249050 161249050 161249050 2 161244738,161246874, 161244895,161249050, 0 LINC01806 unk unk -1,-1, $ zcat hg38.refGene.txt.gz | tail -2 1006 NM_130467 chrX + 55220345 55224108 55220599 55224003 5 55220345,55221374,55221766,55222620,55223986, 55220651,55221463,55221875,55222746,55224108, 0 PAGE5 cmpl cmpl 0,1,0,1,1, 637 NM_001364814 chrY - 6865917 6874027 6866072 6872608 7 6865917,6868036,6868731,6868867,6870005,6872554,6873971, 6866078,6868462,6868776,6868909,6870053,6872620,6874027, 0 AMELY cmpl cmpl 0,0,0,0,0,0,-1,
Where to download reference genome
- UCSC and twoBitToFa to UCSC convert .2bit to fasta.
- hg19 from UCSC (chromosome-wise).
Which human reference genome to use?
http://lh3.github.io/2017/11/13/which-human-reference-genome-to-use (11/13/2017)
GENCODE transcript database
RefSeq categories
See Table 1 of Chapter 18The Reference Sequence (RefSeq) Database.
| Category | Description |
|---|---|
| NC | Complete genomic molecules |
| NG | Incomplete genomic region |
| NM | mRNA |
| NR | ncRNA |
| NP | Protein |
| XM | predicted mRNA model |
| XR | predicted ncRNA model |
| XP | predicted Protein model (eukaryotic sequences) |
| WP | predicted Protein model (prokaryotic sequences) |
UCSC version & NCBI release corresponding
Gene Annotation
- GeneCards
- Genetics Home Reference from National Library of Medicine
- My Cancer Genome
- Cosmic
- annotables R package.
- ensembldb R package
- For gene symbols, there are NCBI and HUGO. See an example from GSE6532 (annot.tam object).
How many DNA strands are there in humans?
- http://www.numberof.net/number-of-dna-strands/
- http://www.answers.com/Q/How_many_DNA_strands_are_there_in_humans
How many base pairs in human
- 3 billion base pairs. https://en.wikipedia.org/wiki/Human_genome
- chromosome 22 has the smallest number of bps (~50 million).
- chromosome 1 has the largest number of bps (245 million base pairs).
- Illumina iGenome Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa file is 3.0GB (so is other genome.fa from human).
Gene, Transcript, Coding/Non-coding exon
- https://hslnews.wordpress.com/2015/07/02/bioinformatics-bite-how-to-find-the-transcription-start-site-of-a-gene/
- According to https://en.wikipedia.org/wiki/Exon, in the human genome only
- 1.1% of the genome is spanned by exons,
- 24% is in introns,
- 75% of the genome being intergenic DNA.
- https://twitter.com/ensembl/status/1276469013707665410?s=20
SNP
Types of SNPs and number of SNPs in each chromosomes
NGS technology
DNA methylation
- GC content = (G+C)/(A+T+G+C) x 100%
- How many CpGs (C follows by G)?
- Analyzing DNA methylation data (part of the book Biomedical Data Science) and the PH525x: Data Analysis for Genomics (edX course). The Github website is on https://github.com/genomicsclass/labs. The source code may not be correct. See also http://www.biostat.jhsph.edu/~iruczins/teaching/kogo/html/ml/week8/methylation.Rmd. The paper Bump hunting to identify differentially methylated regions in epigenetic epidemiology studies the tutorial has mentioned.
devtools::install_github("coloncancermeth","genomicsclass")
library(coloncancermeth) # 485512 x 26
data(coloncancermeth) # load meth (methylation data), pd (sample info ) and gr objects
dim(meth)
dim(pd)
length(gr)
colnames(pd)
table(pd$Status) # 9 normals, 17 cancers
normalIndex <- which(pd$Status=="normal")
cancerlIndex <- which(pd$Status=="cancer")
i=normalIndex[1]
plot(density(meth[,i],from=0,to=1),main="",ylim=c(0,3),type="n")
for(i in normalIndex){
lines(density(meth[,i],from=0,to=1),col=1)
}
### Add the cancer samples
for(i in cancerlIndex){
lines(density(meth[,i],from=0,to=1),col=2)
}
# finding regions of the genome that are different between cancer and normal samples
library(limma)
X<-model.matrix(~pd$Status)
fit<-lmFit(meth,X)
eb <- ebayes(fit)
# plot of the region surrounding the top hit
library(GenomicRanges)
i <- which.min(eb$p.value[,2])
middle <- gr[i,]
Index<-gr%over%(middle+10000)
cols=ifelse(pd$Status=="normal",1,2)
chr=as.factor(seqnames(gr))
pos=start(gr)
plot(pos[Index],fit$coef[Index,2],type="b",xlab="genomic location",ylab="difference")
matplot(pos[Index],meth[Index,],col=cols,xlab="genomic location")
# http://www.ncbi.nlm.nih.gov/pubmed/22422453
# within each chromosome we usually have big gaps creating subgroups of regions to be analyzed
chr1Index <- which(chr=="chr1")
hist(log10(diff(pos[chr1Index])),main="",xlab="log 10 method")
library(bumphunter)
cl=clusterMaker(chr,pos,maxGap=500)
table(table(cl)) ##shows the number of regions with 1,2,3, ... points in them
＃consider two example regions＃
...
Whole Genome Sequencing, Whole Exome Sequencing, Transcriptome (RNA) Sequencing
- http://www.rna-seqblog.com/exome-sequencing-vs-rna-seq-to-identify-coding-region-variants/
- http://www.rna-seqblog.com/combined-use-of-exome-and-transcriptome-sequencing/
- A comparison of whole genome and whole transcriptome sequencing
Sequence + Expression
Integrate RNA-Seq and DNA-Seq
- Integrated RNA and DNA sequencing reveals early drivers of metastatic breast cancer by Perou. An R code is provided.
Integrate/combine Omics
- OmicsPLS & the paper in BMC Bioinformatics 2018
- MultiAssayExperiment & the paper in AACR 2017
- mixOmics, http://mixomics.org/, the paper on PLOS 2017
- The impact of different sources of heterogeneity on loss of accuracy from genomic prediction models Biostatistics 2018
Gene expression
Expression level is the amount of RNA in cell that was transcribed from that gene. Slides from Alyssa Frazee.
Quantile normalization
- https://en.wikipedia.org/wiki/Quantile_normalization
- Question: Quantile normalizing prior to or after TPM scaling?
- When to use Quantile Normalization? and its R package quantro
- normalize.quantiles() from preprocessCore package. Note for ties, the average is used in normalize.quantiles(), ((4.666667 + 5.666667) / 2) = 5.166667.
source('http://bioconductor.org/biocLite.R') biocLite('preprocessCore') #load package library(preprocessCore) #the function expects a matrix #create a matrix using the same example mat <- matrix(c(5,2,3,4,4,1,4,2,3,4,6,8), ncol=3) mat # [,1] [,2] [,3] #[1,] 5 4 3 #[2,] 2 1 4 #[3,] 3 4 6 #[4,] 4 2 8 #quantile normalisation normalize.quantiles(mat) # [,1] [,2] [,3] #[1,] 5.666667 5.166667 2.000000 #[2,] 2.000000 2.000000 3.000000 #[3,] 3.000000 5.166667 4.666667 #[4,] 4.666667 3.000000 5.666667
Merging two gene expression studies; batch effect
- Statistics for Genomic Data Science (Coursera) and https://github.com/jtleek/genstats
- Alternative empirical Bayes models for adjusting for batch effects in genomic studies Zhang et al. BMC Bioinformatics 2018. The R package is BatchQC from Bioconductor.
- Combat() function in sva package from Bioconductor.
- It can remove both known batch effects and other potential latent sources of variation.
- The tutorial includes information on (1) how to estimate the number of latent sources of variation, (2) how to apply the sva package to estimate latent variables such as batch effects, (3) how to directly remove known batch effects using the ComBat function, (4) how to perform differential expression analysis using surrogate variables either directly or with thelimma package, and (4) how to apply “frozen” sva to improve prediction and clustering.
- Figure S1 shows the principal component analysis (PCA) before and after batch effect correction for training and validation datasets
- Tutorial example to remove the batch effect
library(sva) library(bladderbatch) data(bladderdata) pheno = pData(bladderEset) edata = exprs(bladderEset) batch = pheno$batch modcombat = model.matrix(~1, data=pheno) combat_edata = ComBat(dat=edata, batch=batch, mod=modcombat, par.prior=TRUE, prior.plots=FALSE) # This returns an expression matrix, with the same dimensions # as your original dataset. # By default, it performs parametric empirical Bayesian adjustments. # If you would like to use nonparametric empirical Bayesian adjustments, # use the par.prior=FALSE option (this will take longer).
- Merging two gene-expression studies via cross-platform normalization by Shabalin et al, Bioinformatics 2008. This method (called Cross-Platform Normalization/XPN)was used by Ternès Biometrical Journal 2017.
- Batch effect removal methods for microarray gene expression data integration: a survey by Lazar et al, Bioinformatics 2012. The R package is inSilicoMerging which has been removed from Bioconductor 3.4.
- Question: Combine hgu133a&b and hgu133plus2. Adjusting batch effects in microarray expression data using empirical Bayes methods
- removeBatchEffect() from limma package
- Batch effects and GC content of NGS by Michael Love
Fusion gene
https://en.wikipedia.org/wiki/Fusion_gene
Structural variation
- https://en.wikipedia.org/wiki/Structural_variation
- https://www.ncbi.nlm.nih.gov/dbvar/content/overview/
- Detection of complex structural variation from paired-end sequencing data Joseph G. Arthur, 2018
LUMPY, DELLY, ForestSV, Pindel, breakdancer , SVDetect.
RNASeq + ChipSeq
Labs
Biowulf2 at NIH
- Main site: http://hpc.nih.gov
- User guide: https://hpc.nih.gov/docs/user_guides.html
- Unlock account (60 days inactive) https://hpc.nih.gov/dashboard/
- Transitioning from PBS to Slurm: https://hpc.nih.gov/docs/pbs2slurm.html
- Job Submission 'cheat sheet': https://hpc.nih.gov/docs/biowulf2-handout.pdf
- STAR: https://hpc.nih.gov/apps/STAR.html
BamHash
Hash BAM and FASTQ files to verify data integrity. The C++ code is based on OpenSSL and seqan libraries.
Reproducibility
- Improving reproducibility in computational biology research
- SnakeChunks: modular blocks to build Snakemake workflows for reproducible NGS analyses by Claire Rioualen et al, 2017.
Selected Papers
- Testing for association between RNA-Seq and high-dimensional data and the Bioconductor package globalSeq.
- The FDA’s Experience with Emerging Genomics Technologies—Past, Present, and Future
- Differential analysis of gene regulation at transcript resolution with RNA-seq Trapnell et al, Nature Biotechnology 31, 46–53 (2013)
- A Study of TP53 RNA Splicing Illustrates Pitfalls of RNA-seq Methodology
- Top RNA-Seq Articles – 2016 from RNA-Seq blog
- Multivariate association analysis with somatic mutation data by He 2017 Biometrics.
- RASflow – An RNA-Seq Analysis Workflow with Snakemake Zhang 2019
- A Survey of Bioinformatics Database and Software Usage through Mining the Literature
Pictures
https://www.flickr.com/photos/genomegov
用DNA做身分鑑識
如何自学入门生物信息学
https://zhuanlan.zhihu.com/p/32065916
Staying current
Staying Current in Bioinformatics & Genomics: 2017 Edition
Papers
Common issues in algorithmic bioinformatics papers
What are some common issues I find when reviewing algorithmic bioinformatics conference papers?
Precision Medicine courses
Personalized medicine
- F.D.A. Approves First Gene-Altering Leukemia Treatment
- The FDA Just Approved a New Way of Fighting (lymphoma) Cancer Using Personalized Gene Therapy
- What is the difference between precision medicine and personalized medicine? What about pharmacogenomics? "personalized medicine" is an older term. Help Me Understand Genetics
Cancer and gene markers
- Colorectal cancer patients without KRAS mutations have far better outcomes with EGFR treatment than those with KRAS mutations.
- Two EGFR inhibitors, cetuximab and panitumumab are not recommended for the treatment of colorectal cancer in patients with KRAS mutations in codon 12 and 13.
- Breast cancer.
The shocking truth about space travel
bioSyntax: syntax highlighting for computational biology
https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-018-2315-y
Deep learning
Deep learning: new computational modelling techniques for genomics
Terms
基因结构
https://zhuanlan.zhihu.com/p/49601643
Epidemiology
Epidemiology for the uninitiated
in vitro
- https://en.wikipedia.org/wiki/In_vitro
- https://en.wikipedia.org/wiki/In_silico
- https://en.wikipedia.org/wiki/RNA-Seq
RNA sequencing 101
Web
- Introduction to Biology - The Secret of Life from edX. It's one of the top 100 best free online courses posted by ClassCentral.
- Introduction to RNA-Seq including biology overview (DNA, Alternative splcing, mRNA structure, human genome) and sequencing technology.
- RNA Sequencing 101 by Agilent Technologies. Includes the definition of sequencing depth (number of reads per sample) and coverage (number of reads/locus).
- Where do we get reads(A,C,T,G) from sample RNA? See page 12 of this pdf from Colin Dewey in U. Wisc.
- Quantification of RNA-Seq data (see the above pdf)
- Convert read counts into expression: RPKM (see the above pdf)
- RPKM and FPKM (Data analysis of RNA-seq from new generation sequencing by 張庭毓) RNA-Seq資料分析研討會與實作課程 / RNA Seq定序 / 次世代定序(NGS) / 高通量基因定序 分析.
- First vs Second generation sequence.
- The Simple Fool’s Guide to Population Genomics via RNASeq: An Introduction to HighThroughput Sequencing Data Analysis. This covers QC, De novo assembly, BLAST, mapping reads to reference sequences, gene expression analysis and variant (SNP) detection.
- An Introduction to Bioinformatics Resources and their Practical Applications from NIH library Bioinformatics Support Program.
- Teaching material from rnaseqforthenextgeneration.org which includes Designing RNA-Seq experiments, Processing RNA-Seq data, and Downstream analyses with RNA-Seq data.
Books
- RNA-seq Data Analysis: A Practical Approach. The pdf version is available on slideshare.net.
- Statistical Analysis of Next Generation Sequencing Data
- Modern Statistics for Modern Biology (free, see Statistics books). Plots are based on ggplot2.
- Learning Bioinformatics At Home - some resources gathered by the Harvard Informatics group.
- Install all required packages using the R script gives me errors. It would be great if there is a Docker image. Another solution is to run install.packages(pkgsToInstall, Ncpus=4) manually.
strand-specific vs non-strand specific experiment
- http://seqanswers.com/forums/showthread.php?t=28025. According to this message and the article (under the paragraph of Read counting) from PLOS, most of the RNA-seq protocols that are used nowadays are not strand-specific.
- http://biology.stackexchange.com/questions/1958/difference-between-strand-specific-and-not-strand-specific-rna-seq-data
- https://www.biostars.org/p/61625/ how to find if this public rnaseq data are prepared by strand-specific assay?
- https://www.biostars.org/p/62747/ Discussion of using IGV to view strand-specific coverage. See also the similar posts on the right hand side.
- https://www.biostars.org/p/44319/ How To Find Stranded Rna-Seq Experiments Data. The text dUTP 2nd strand marking includes a link to stranded rna-seq data.
- forward (+)/ reverse(-) strand in GAlignments objects (p68 of the pdf manual and page 7 of sam format specification.
Understand this info is necessary when we want to use summarizeOverlaps() function (GenomicAlignments) or htseq-count python program to get count data.
This post mentioned to use infer_experiment.py script to check whether the rna-seq run is stranded or not.
The rna-seq experiment used in this tutorial is not stranded-specific.
FASTQ
- FASTQ=FASTA + Qual. FASTQ format is a text-based format for storing both a biological sequence (usually nucleotide sequence) and its corresponding quality scores.
Phred quality score
q = -10log10(p) where p = error probability for the base.
| q | error probability | base call accuracy |
|---|---|---|
| 10 | 0.1 | 90% |
| 13 | 0.05 | 95% |
| 20 | 0.01 | 99% |
| 30 | 0.001 | 99.9% |
| 40 | 0.0001 | 99.99% |
| 50 | 0.00001 | 99.999% |
FASTA
fasta/fa files can be used as reference genome in IGV. But we cannot load these files in order to view them.
Download sequence files
- ftp://ftp.ncbi.nih.gov/genomes/H_sapiens/RNA/ Assembled genome sequence and annotation data for RefSeq genome assemblies
- ftp://ftp.ncbi.nlm.nih.gov/refseq/H_sapiens/
Compute the sequence length of a FASTA file
https://stackoverflow.com/questions/23992646/sequence-length-of-fasta-file
awk '/^>/ {if (seqlen){print seqlen}; print ;seqlen=0;next; } { seqlen += length($0)}END{print seqlen}' file.fa
head -2 file.fa | \
awk '/^>/ {if (seqlen){print seqlen}; print ;seqlen=0;next; } { seqlen += length($0)}END{print seqlen}' | \
tail -1
FASTA <=> FASTQ conversion
According to this post,
- FastA are text files containing multiple DNA* seqs each with some text, some part of the text might be a name.
- FastQ files are like fasta, but they also have quality scores for each base of each seq, making them appropriate for reads from an Illumina machine (or other brands)
Convert FASTA to FASTQ without quality scores
Biostars. For example, the bioawk by lh3 (Heng Li) worked.
Convert FASTA to FASTQ with quality score file
See the links on the above post.
Convert FASTQ to FASTA using Seqtk
Use the Seqtk program; see this post.
The Seqtk program by lh3 can be used to sample reads from a fastq file including paired-end; see this post.
RPKM (Mortazavi et al. 2008) and cpm (counts per million)
Reads per Kilobase of Exon per Million of Mapped reads.
- rpkm function in edgeR package.
- RPKM function in easyRNASeq package.
- TMM > cpm > log2 transformation on the paper Gene expression profiling of 1200 pancreatic ductal adenocarcinoma reveals novel subtypes
- Gene expression quantification from RNA-Seq wikipedia page
- Sequencing depth/coverage: the total number of reads generated in a single experiment is typically normalized by converting counts to fragments, reads, or counts per million mapped reads (FPM, RPM, or CPM).
- Gene length: FPKM, TPM. Longer genes will have more fragments/reads/counts than shorter genes if transcript expression is the same. This is adjusted by dividing the FPM by the length of a gene, resulting in the metric fragments per kilobase of transcript per million mapped reads (FPKM). When looking at groups of genes across samples, FPKM is converted to transcripts per million (TPM) by dividing each FPKM by the sum of FPKMs within a sample.
- Difference between CPM and TPM and which one for downstream analysis?. CPM is basically depth-normalized counts whereas TPM is length normalized (and then normalized by the length-normalized values of the other genes).
- The RNA-seq abundance zoo. The counts per million (CPM) metric takes the raw (or estimated) counts, and performs the first type of normalization I mention in the previous section. That is, it normalized the count by the library size, and then multiplies it by a million (to avoid scary small numbers).
Idea
- The more we sequence, the more reads we expect from each gene. This is the most relevant correction of this method.
- Longer transcript are expected to generate more reads. The latter is only relevant for comparisons among different genes which we rarely perform!. As such, the DESeq2 only creates a size factor for each library and normalize the counts by dividing counts by a size factor (scalar) for each library. Note that: H0: mu1=mu2 is equivalent to H0: c*mu1=c*mu2 where c is gene length.
Calculation
- Count up the total reads in a sample and divide that number by 1,000,000 – this is our “per million” scaling factor.
- Divide the read counts by the “per million” scaling factor. This normalizes for sequencing depth, giving you reads per million (RPM)
- Divide the RPM values by the length of the gene, in kilobases. This gives you RPKM.
Formula
RPKM = (10^9 * C)/(N * L), with C = Number of reads mapped to a gene N = Total mapped reads in the experiment L = gene length in base-pairs for a gene
source("http://www.bioconductor.org/biocLite.R")
biocLite("edgeR")
library(edgeR)
set.seed(1234)
y <- matrix(rnbinom(20,size=1,mu=10),5,4)
[,1] [,2] [,3] [,4]
[1,] 0 0 5 0
[2,] 6 2 7 3
[3,] 5 13 7 2
[4,] 3 3 9 11
[5,] 1 2 1 15
d <- DGEList(counts=y, lib.size=1001:1004)
# Note that lib.size is optional
# By default, lib.size = colSums(counts)
cpm(d) # counts per million
Sample1 Sample2 Sample3 Sample4
1 0.000 0.000 4985.045 0.000
2 5994.006 1996.008 6979.063 2988.048
3 4995.005 12974.052 6979.063 1992.032
4 2997.003 2994.012 8973.081 10956.175
5 999.001 1996.008 997.009 14940.239
> cpm(d,log=TRUE)
Sample1 Sample2 Sample3 Sample4
1 7.961463 7.961463 12.35309 7.961463
2 12.607393 11.132027 12.81875 11.659911
3 12.355838 13.690089 12.81875 11.129470
4 11.663897 11.662567 13.17022 13.451207
5 10.285119 11.132027 10.28282 13.890078
d$genes$Length <- c(1000,2000,500,1500,3000)
rpkm(d)
Sample1 Sample2 Sample3 Sample4
1 0.0000 0.000 4985.0449 0.000
2 2997.0030 998.004 3489.5314 1494.024
3 9990.0100 25948.104 13958.1256 3984.064
4 1998.0020 1996.008 5982.0538 7304.117
5 333.0003 665.336 332.3363 4980.080
> cpm
function (x, ...)
UseMethod("cpm")
<environment: namespace:edgeR>
> showMethods("cpm")
Function "cpm":
<not an S4 generic function>
> cpm.default
function (x, lib.size = NULL, log = FALSE, prior.count = 0.25,
...)
{
x <- as.matrix(x)
if (is.null(lib.size))
lib.size <- colSums(x)
if (log) {
prior.count.scaled <- lib.size/mean(lib.size) * prior.count
lib.size <- lib.size + 2 * prior.count.scaled
}
lib.size <- 1e-06 * lib.size
if (log)
log2(t((t(x) + prior.count.scaled)/lib.size))
else t(t(x)/lib.size)
}
<environment: namespace:edgeR>
> rpkm.default
function (x, gene.length, lib.size = NULL, log = FALSE, prior.count = 0.25,
...)
{
y <- cpm.default(x = x, lib.size = lib.size, log = log, prior.count = prior.count)
gene.length.kb <- gene.length/1000
if (log)
y - log2(gene.length.kb)
else y/gene.length.kb
}
<environment: namespace:edgeR>
Here for example the 1st sample and the 2nd gene, its rpkm value is calculated as
# step 1: 6/(1.0e-6 *1001) = 5994.006 # cpm, compute column-wise # step 2: 5994.006/ (2000/1.0e3) = 2997.003 # rpkm, compute row-wise # Another way # step 1 (RPK) 6/ (2000/1.0e3) = 3 # step 2 (RPKM) 3/ (1.0e-6 * 1001) = 2997.003
Critics
Consider the following example: in two libraries, each with one million reads, gene X may have 10 reads for treatment A and 5 reads for treatment B, while it is 100x as many after sequencing 100 millions reads from each library. In the latter case we can be much more confident that there is a true difference between the two treatments than in the first one. However, the RPKM values would be the same for both scenarios. Thus, RPKM/FPKM are useful for reporting expression values, but not for statistical testing!
(another critic) Union Exon Based Approach
http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0141910
In general, the methods for gene quantification can be largely divided into two categories: transcript-based approach and ‘union exon’-based approach.
It was found that the gene expression levels are significantly underestimated by ‘union exon’-based approach, and the average of RPKM from ‘union exons’-based method is less than 50% of the mean expression obtained from transcript-based approach.
FPKM (Trapnell et al. 2010)
Fragment per Kilobase of exon per Million of Mapped fragments (Cufflinks). FPKM is very similar to RPKM. RPKM was made for single-end RNA-seq, where every read corresponded to a single fragment that was sequenced. FPKM was made for paired-end RNA-seq. With paired-end RNA-seq, two reads can correspond to a single fragment, or, if one read in the pair did not map, one read can correspond to a single fragment. The only difference between RPKM and FPKM is that FPKM takes into account that two reads can map to one fragment (and so it doesn’t count this fragment twice).
RPKM, FPKM, TPM and DESeq
- The youtube video is on here. TPM 比 RPKM/FPKM 好因為 total reads in each experiments are the same.
- Relationship (formula) between RPKM and TPM
- Differences
- The main difference between RPKM and FPKM is that the former is a unit based on single-end reads, while the latter is based on paired-end reads and counts the two reads from the same RNA fragment as one instead of two.
- The difference between RPKM/FPKM and TPM is that the former calculates sample-scaling factors before dividing read counts by gene lengths, while the latter divides read counts by gene lengths first and calculates samples calling factors based on the length-normalized read counts.
- If researchers would like to interpret gene expression levels as the proportions of RNA molecules from different genes in a sample,
- 为什么说FPKM和RPKM都错了？
- Measurement of mRNA abundance using RNA-seq data: RPKM measure is inconsistent among samples (TPM method) by Wagner 2012.
- Between samples normalization.
- How to Normalize Salmon TPM output?
- How do I normalize for my RNA-seq data across different samples in different conditions. Using DESeq2 (paper, estimateSizeFactors()). DESeq paper by Anders 2010 where the NB model and size factor was first used.
> dds <- makeExampleDESeqDataSet(m=4) > head(counts(dds)) sample1 sample2 sample3 sample4 gene1 0 1 1 6 gene2 2 0 0 3 gene3 18 9 19 12 gene4 12 25 13 13 gene5 22 26 10 8 gene6 6 5 8 6 > dds <- estimateSizeFactors(dds) > head(counts(dds)) sample1 sample2 sample3 sample4 gene1 0 1 1 6 gene2 2 0 0 3 gene3 18 9 19 12 gene4 12 25 13 13 gene5 22 26 10 8 gene6 6 5 8 6 > head(counts(dds, normalized=TRUE)) sample1 sample2 sample3 sample4 gene1 0.00000 0.9654796 0.9858756 5.732657 gene2 1.96066 0.0000000 0.0000000 2.866328 gene3 17.64594 8.6893164 18.7316365 11.465314 gene4 11.76396 24.1369899 12.8163829 12.420756 gene5 21.56726 25.1024695 9.8587560 7.643542 gene6 5.88198 4.8273980 7.8870048 5.732657
- TPM (generated by RSEM or Salmon or Kallisto) has been suggested as a better unit than RPKM/FPKM. But it cannot be used to do comparison between samples.
- TPM and FPKM from array suite wiki
- What the FPKM? A review of RNA-Seq expression units. R code for computing effective counts, TPM, and FPKM/
- In RNA-Seq, 2 != 2: Between-sample normalization. Among the most popular and well-accepted BSN (between-sample normalization) methods are TMM and DESeq normalization.
P -- per K -- kilobase (related to gene length) M -- million (related to sequencing depth)
- How to convert transcript level TPM to gene level TPM ?
- Calculating TPM from featureCounts output
- Comparison of common normalization methods from a workshop from Harvard Chan Bioinformatics Core.
TMM (Robinson and Oshlack, 2010)
Trimmed Means of M values (EdgeR).
Sample size
Coverage
- Recommended Coverage and Read Depth for NGS Applications from @Genohub.
- bedtools. The bedtools is now hosted on github
- https://github.com/alyssafrazee/polyester
~20x coverage ----> reads per transcript = transcriptlength/readlength * 20
- Page 18 of this RNA-Seq 101 from Agilent or Estimating Sequencing Coverage from Illumina.
C = L N / G
where L=read length, N =number of reads and G=haploid genome length. So, if we take one lane of single read human sequence with v3 chemistry, we get C = (100 bp)*(189×10^6)/(3×10^9 bp) = 6.3. This tells us that each base in the genome will be sequenced between six and seven times on average.
- coverage() function in IRanges package.
- bamcov - Quickly calculate and visualize sequence coverage in alignment files
- Coverage and read depth
- What Is The Sequencing 'Depth' ? Coverage = (total number of bases generated) / (size of genome sequenced). So a 30x coverage means, on an average each base has been read by 30 sequences.
- Sequencing depth and coverage: key considerations in genomic analyses
- Qualimap
- https://www.biostars.org/p/5165/
# Assume the bam file is sorted by chromosome location
# took 40 min on 5.8G bam file. samtools depth has no threads option:(
# it is not right since it only account for regions that were covered with reads
samtools depth *bamfile* | awk '{sum+=$3} END { print "Average = ",sum/NR}' # maybe 42
# The following is the right way! The result matches with Qualimap program.
samtools depth -a *bamfile* | awk '{sum+=$3} END { print "Average = ",sum/NR}' # maybe 8
# OR
LEN=`samtools view -H bamfile | grep -P '^@SQ' | cut -f 3 -d ':' | awk '{sum+=$1} END {print sum}'` # 3095693981
SUM=`samtools depth bamfile | awk '{sum+=$3} END { print "Sum = ", sum}'` # 24473867730
echo $(( $LEN/$SUM ))
SAM/Sequence Alignment Format and BAM format specification
- https://samtools.github.io/hts-specs/SAMv1.pdf and samtools webpage.
- http://genome.sph.umich.edu/wiki/SAM
Single-end, pair-end, fragment, insert size
Germline vs Somatic mutation
Germline: inherit from parents. See the Wikipedia page.
Driver vs passenger mutation
https://en.wikipedia.org/wiki/Somatic_evolution_in_cancer
Nonsynonymous mutation
It is related to the genetic code, Wikipedia. There are 20 amino acids though there are 64 codes.
See
- http://evolution.about.com/od/Overview/a/Synonymous-Vs-Nonsynonymous-Mutations.htm
- https://en.wikipedia.org/wiki/Nonsynonymous_substitution
- Understanding SNPs and INDELs in microbial genomes
- An example from https://en.wikipedia.org/wiki/Silent_mutation
- nonsynonymous: ATG to GTG mutation (AUG = Met, GUG = Val)
- synonymous: CAT to CAC mutation (CAU = His, CAC = His)
isma: analysis of mutations detected by multiple pipelines
isma: an R package for the integrative analysis of mutations detected by multiple pipelines
Tumor mutational burden
- Significance of tumor mutation burden combined with immune infiltrates in the progression and prognosis of ovarian cancer
- maftools Summarize, Analyze and Visualize MAF Files
- 2449 result(s) for 'TMB' from BMC
Missense variants
aminoacid changing variants
Alternative and differential splicing
Best practices and appropriate workflows to analyse alternative and differential splicing
Allele vs Gene
http://www.diffen.com/difference/Allele_vs_Gene
- A gene is a stretch of DNA or RNA that determines a certain trait.
- Genes mutate and can take two or more alternative forms; an allele is one of these forms of a gene. For example, the gene for eye color has several variations (alleles) such as an allele for blue eye color or an allele for brown eyes.
- An allele is found at a fixed spot on a chromosome?
- Chromosomes occur in pairs so organisms have two alleles for each gene — one allele in each chromosome in the pair. Since each chromosome in the pair comes from a different parent, organisms inherit one allele from each parent for each gene. The two alleles inherited from parents may be same (homozygous) or different (heterozygotes).
Locus
https://en.wikipedia.org/wiki/Locus_(genetics)
Haplotypes
- http://www.brown.edu/Research/Istrail_Lab/proj_cmsh.php
- http://www.nature.com/nri/journal/v5/n1/fig_tab/nri1532_F2.html
- https://www.sciencenews.org/article/seeking-genetic-fate
- http://www.medscape.com/viewarticle/553400_3
Base quality, Mapping quality, Variant quality
- Fastq base quality: https://en.wikipedia.org/wiki/FASTQ_format
- Mapping quality: http://genome.cshlp.org/content/18/11/1851.long
- Variant quality: http://www.ncbi.nlm.nih.gov/pubmed/21903627
Mapping quality (MAPQ) vs Alignment score (AS)
http://seqanswers.com/forums/showthread.php?t=66634 & SAM format specification
- MAPQ (5th column): MAPping Quality. It equals −10 log10 Pr{mapping position is wrong} (defined by SAM documentation), rounded to the nearest integer. A value 255 indicates that the mapping quality is not available. MAPQ is a metric that tells you how confident you can be that the read comes from the reported position. So given 1000 reads, for example, read alignments with mapping quality being 30, one of them will be wrong in average (10^(30/-10)=.001). Another example, if MAPQ=70, then the probability mapping position is wrong is 10^(70/-10)=1e-7. We can use 'samtools view -q 30 input.bam' to keep reads with MAPQ at least 30. Users should refer to the alignment program for the 'MAPQ' value it uses.
- AS (optional, 14th column in my case): Alignment score is a metric that tells you how similar the read is to the reference. AS increases with the number of matches and decreases with the number of mismatches and gaps (rewards and penalties for matches and mismatches depend on the scoring matrix you use)
Note:
- MAPQ scores produced by the aligners typically involves the alignment score and other information.
- You can have high AS and low MAPQ if the read aligns perfectly at multiple positions, and you can have low AS and high MAPQ if the read aligns with mismatches but still the reported position is still much more probable than any other.
- You probably want to filter for MAPQ, but "good" alignment may refer to AS if what you care is similarity between read and reference.
- MAPQ values are really useful but their implementation is a mess by Simon Andrews
gene's isoform
- Bioinformatics Tools for RNA-seq Gene and Isoform Quantification 2016. Figure 1 gives an illustration how isoform can affect the computation of RPKM.
- Evaluation and comparison of computational tools for RNA-seq isoform quantification BMC Genomics 2017.
- Differential Isoform Expression With RNA-Seq: Are We Really There Yet?
- An example of unidentified transcript-level estimates while gene-level estimation is still possible. Figure 1C in Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences
Other software
Partek
- Partek Flow Software http://youtu.be/-6aeQPOYuHY
dCHIP
MeV
MeV v4.8 (11/18/2011) allows annotation from Bioconductor
IPA from Ingenuity
Login: There are web started version https://analysis.ingenuity.com/pa and Java applet version https://analysis.ingenuity.com/pa/login/choice.jsp. We can double click the file <IpaApplication.jnlp> in my machine's download folder.
Features:
- easily search the scientific literature/integrate diverse biological information.
- build dynamic pathway models
- quickly analyze experimental data/Functional discovery: assign function to genes
- share research and collaborate. On the other hand, IPA is web based, so it takes time for running analyses. Once submitted analyses are done, an email will be sent to the user.
Start Here
Expression data -> New core analysis -> Functions/Diseases -> Network analysis
Canonical pathways |
| |
Simple or advanced search --------------------+ |
| |
v |
My pathways, Lists <------+
^
|
Creating a custom pathway --------------------+
Resource:
- http://bioinformatics.mdanderson.org/MicroarrayCourse/Lectures09/Pathway%20Analysis.pdf
- http://libguides.mit.edu/content.php?pid=14149&sid=843471
- http://people.mbi.ohio-state.edu/baguda/PathwayAnalysis/
- IPA 5.5 manual http://people.mbi.ohio-state.edu/baguda/PathwayAnalysis/ipa_help_manual_5.5_v1.pdf
- Help and supports
- Tutorials which includes
- Search for genes
- Analysis results
- Upload and analyze example data
- Upload and analyze your own expression data
- Visualize connections among genes
- Learn more special features
- Human isoform view
- Transcription factor analysis
- Downstream effects analysis
Notes:
- The input data file can be an Excel file with at least one gene ID and expression value at the end of columns (just what BRB-ArrayTools requires in general format importer).
- The data to be uploaded (because IPA is web-based; the projects/analyses will not be saved locally) can be in different forms. See http://ingenuity.force.com/ipa/articles/Feature_Description/Data-Upload-definitions. It uses the term Single/Multiple Observation. An Observation is a list of molecule identifiers and their corresponding expression values for a given experimental treatment. A dataset file may contain a single observation or multiple observations. A Single Observation dataset contains only one experimental condition (i.e. wild-type). A Multiple Observation dataset contains more than one experimental condition (i.e. a time course experiment, a dose response experiment, etc) and can be uploaded into IPA in a single file (e.g. Excel). A maximum of 20 observations in a single file may be uploaded into IPA.
- The instruction http://ingenuity.force.com/ipa/articles/Feature_Description/Data-Upload-definitions shows what kind of gene identifier types IPA accepts.
- In this prostate example data tutorial, the term 'fold change' was used to replace log2 gene expression. The tutorial also uses 1.5 as the fold change expression cutoff.
- The gene table given on the analysis output contains columns 'Fold change', 'ID', 'Notes', 'Symbol' (with tooltip), 'Entrez Gene Name', 'Location', 'Types', 'Drugs'. See a screenshot below.
Screenshots:
File:IngenuityAnalysisOutput.png
DAVID Bioinformatics Resource
It offers an integrated annotation combining gene ontology, pathways and protein annotations.
It can be used to identify the pathways associated with a set of genes; e.g. this paper.
GOTrapper
GOTrapper: a tool to navigate through branches of gene ontology hierarchy
qpcR
Model fitting, optimal model selection and calculation of various features that are essential in the analysis of quantitative real-time polymerase chain reaction (qPCR).
GSEA
- http://www.broadinstitute.org/gsea/doc/desktop_tutorial.jsp
- http://www.hmwu.idv.tw/web/CourseSMDA/MADA/Hank_MicroarrayDataAnalysis-GSEA-20110616.pdf
- MGSEA– a multivariate Gene set enrichment analysis