Wednesday, February 15, 2017

gene expression profile outlier analysis


1) Packages that makes parametric assumption - mostly for microarray gene expression data.

MOST, ZOPET, OPPAR, mCOPA, PPST, OS, ORT.

2) non- parametric

coGPS. From the first glance look very useful.

#1 - non parametric
#2 - allows integration of multiple data types like expression, methylation and copy number.
#3 - permutation based p-value (based on D. Ghosh's work)
#4 -  Both gene and patience wise outlier detection.

Wednesday, November 11, 2015

Installing MACS2 without hiccups

Installing MACS2 by Tao Liu

I was interested in identifying differential peak from ChIP-Seq data without replicates.
So, as a happy user of MACS, I wanted to use MACS2.

Tips:

1) Install numpy
2) Install Cython (this would remove Prob.c error)
3) In "setup_w_cython.py" script - change -Ofast to -O3 to make sure GCC is not offended.

It should compile smooth.

Here is how to run MACS2
https://github.com/taoliu/MACS/wiki/Call-differential-binding-events

Will update on result quality.

 

Friday, December 20, 2013

Turning R string in to a variable

http://cran.r-project.org/doc/FAQ/R-FAQ.html#How-can-I-turn-a-string-into-a-variable_003f

"get" and "assign" allows "paste" to generate variables. 

Infinitely Abusable !

Saturday, November 30, 2013

Setting up High throughput RNA-Seq analysis pipeline

This email provides general discussion of the workflow and tools for analysis of RNA-Seq data. This are not the only tools but things that work.

A general requirement is install of recent version of Java, R and Bioconductor packages.

1) Fastq file(s).

FASTQ format stores sequences and Phred quality scores in a single file.
http://en.wikipedia.org/wiki/FASTQ_format

(As an example) Assuming a Illumina Hi-Seq run of 8 lanes, RNA-Seq data (reads) for up to 32 samples can be generated in 1 run. This results in upto 4 samples in a single lane. After de-multiplexing (if any) based on sequence barcode, one normally receives one or more FASTQ file for a given sample (single experiment).

If the paired-end (PE) sequencing was run the file will have R1 and R2 in the name.

One sample can have more than 1 fastq file and it can be merged with simple PERL script to concatnet them. R1 and R2 files should be merged separately.

2) Alignment to the reference genome assembly. (SAM and BAM files)

The first step in the analysis is to align the sequence reads in the fastq files to the genome.
There are three repositories to download the genome assemblies the UCSC, NCBI and ENSEMBL.

UCSC : http://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/ chromFaMasked.tar.gz provides masked chromosome fasta files. Similarly, ENSEMBL fasta files can be downloaded from the ENSEMBL website (www.ensembl.org)

The UCSC assmebly is called "hg19" and ENSEMBL assembly is called "GRCh37". There is NO difference between them in terms of nucleotide sequences. (annotations differ)

While many tools available for doing the alignment Bowtie (http://bowtie-bio.sourceforge.net/index.shtml), Tophat (http://tophat.cbcb.umd.edu/) and BWA (http://bio-bwa.sourceforge.net/) are most frequently used. You may not need to download genome sequences from UCSC or ENSEMBL website. The websites of the aligners provides assemblies and index files
specifically designed to run with the aligners.

The output of the alignment is SAM files and its binary version is called BAM. see http://samtools.sourceforge.net/

3) QC for the sequencing reads and RNA-Seq specific QC

FASTQC allows for quality control of raw sequencing reads. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/

RNA-SeQC allows for quality control specific of RNA-seq experiments, providing quality metric for transcriptome analysis.
http://www.broadinstitute.org/cancer/cga/rna-seqc

There are other packages available but these two are currently the best.

RNA-SeQC allows marking and removal of duplicate reads. The BAM files from which the duplicate reads are removed is the input of the further analysis.


4) Transcript assembly and Counting of reads for transcriptome.

This is a RNA-Seq analysis specific step. This step is a basic requirement for down stream analysis like differential expression analysis. requires genomic annotation information that is available
from GFF format files. The two places for genomic annotation are UCSC and ENSEMBL. I prefer ENSEMBL annotations but it is nice to use the ENCODE annotations currently used by TCGA.

The softwares are Cufflinks (http://cufflinks.cbcb.umd.edu/), HTSeq
(http://www-huber.embl.de/users/anders/HTSeq/doc/overview.html) and
easyRNASeq (http://www.bioconductor.org/packages/release/bioc/html/easyRNASeq.html)

Various quantitation metrices like RPKM, FPFM (for paired end data), RSEM and raw counts are proposed for RNA-seq data. Current opinion is that RSEM (http://deweylab.biostat.wisc.edu/rsem/) is useful for within sample comparison and raw counts are suitable for differential expression analysis.

5) Normalization, Filtering and Differential Expression

The most commonly used packages are DESeq , DESeq2 and edgeR. There are others but I used mostly DESeq and edgeR for they allow incorporation of experiment design for fitting linear models 
and analysis of variance (ANOVA).

Filtering is important to increase power of the differential expression analysis. See discussion on Independent filtering by the "genefilter" and HTSFilter Bioconductor packages.

6) Visualization

The CummerBund package (http://compbio.mit.edu/cummeRbund/) associated with cufflinks allows for visualization of transcript data. RSEM (http://deweylab.biostat.wisc.edu/rsem/) is useful for within sample comparison and visualization.

7) Functional enrichment and pathway analysis.

No clear winner here. However I use GSVA, GSEA, DAVID and other packages.

Thursday, August 08, 2013

RPKM with easyRNASeq

Seems there is a bug in the RPKM function in easyRNASeq. Trying to generate RPKM values from counts results in the matrix with "NA".

The solution from seqanswers works :
http://seqanswers.com/forums/showthread.php?t=29909

It works for the other version of the package too.

Saturday, June 08, 2013

Nice worked out example on analyzing array data using bioconductor.

http://bioinformatics.knowledgeblog.org/2011/06/20/analysing-microarray-data-in-bioconductor/

Monday, April 15, 2013

R/Bioconductor packages for VennDiagram

gplots::venn
venneuler
limma::vennDiagram
vennDiagram

the paper by Hanbo Chen and Paul Boutros give good comparison of features of multiple tools. 
http://www.biomedcentral.com/1471-2105/12/35

Wednesday, January 02, 2013

Experimental Design for ANOVA

A good introduction is here

graphpad website for ANOVA introduction

The term repeated-measures refers to an experiment that collects multiple measurements from each subject. The analysis of repeated measures data is identical to the analysis of randomized block experiments that use paired or matched subjects. Prism can calculate repeated-measures two-way ANOVA when either one of the factors are repeated or matched (mixed model) or when both factors are.
One data table can correspond to four experimental designs
Prism uses a unique way to enter data. You use rows and columns to designate the different groups (levels) of each factor. Each data set (column) represents a different level of one factor, and each row represents a different level of the other factor. You need to decide which factor is defined by rows, and which by columns. Your choice will not affect the ANOVA results, but the choice is important as it affects the appearance of graphs.
The table above shows example data testing the effects of three doses of drugs in control and treated animals.
These data could have come from four distinct experimental designs.
Not repeated measures
The experiment was done with six animals. Each animal was given one of two treatments at one of three doses. The measurement was then made in duplicate. The value at row 1, column A, Y1 (23) came from the same animal as the value at row 1, column A, Y2 (24). Since the matching is within a treatment group, it is a replicate, not a repeated measure. Analyze these data with ordinary two-way ANOVA, not repeated-measures ANOVA.
Matched values are spread across a rows
The experiment was done with six animals, two for each dose. The control values were measured first in all six animals. Then you applied a treatment to all the animals and made the measurement again. In the table above, the value at row 1, column A, Y1 (23) came from the same animal as the value at row 1, column B, Y1 (28). The matching is by row.
Matched values are stacked into a subcolumn
The experiment was done with four animals. First each animal was exposed to a treatment (or placebo). After measuring the baseline data (dose=zero), you inject the first dose and make the measurement again. Then inject the second dose and measure again. The values in the first Y1 column (23, 34, and 43) were repeated measurements from the same animal. The other three subcolumns came from three other animals. The matching was by column.
Repeated measures in both factors
The experiment was done with two animals. First you measured the baseline (control, zero dose). Then you injected dose 1 and made the next measurement, then dose 2 and measured again. Then you gave the animal the experimental treatment, waited an appropriate period of time, and made the three measurements again. Finally, you repeated the experiment with another animal (Y2). So a single animal provided data from both Y1 subcolumns (23, 34, 43 and 28, 41, 56).
When do you specify which design applies to this experiment?
The example above shows that one grouped data set can represent four different experimental designs. You do not distinguish these designs when creating the data table. The data table doesn't "know" wether or not the data are repeated measures. You should take into account experimental design when choosing how to graph the data. And you must take it into account when performing two-way ANOVA. On the first tab of the two-way ANOVA dialog, you'll designate the experimental design.
Lingo: "Repeated measures" vs. "randomized block" experiments
The term repeated measures is appropriate when you made repeated measurements from each subject.
Some experiments involve matching but not repeated measurements. The term randomized-block describes these kinds of experiments. For example, imagine that the three rows were three different cell lines. All the Y1 data came from one experiment, and all the Y2 data came from another experiment performed a month later. The value at row 1, column A, Y1 (23) and the value at row 1, column B, Y1 (28) came from the same experiment (same cell passage, same reagents). The matching is by row.
Randomized block data are analyzed identically to repeated-measures data. Prism always uses the term repeated measures, so you should choose repeated measures analyses when your experiment follows a randomized block design.

Wednesday, September 19, 2012

Adding data points to boxplots with stripchart

pre<-55+rnorm(20)
post<-pre+0.7+rnorm(20)
 
s<-seq(length(pre))
 
par(bty="l")
 
boxplot(pre,post,main="Rawdata",xlab="Time",ylab="Measure",names=c("pre","post"), 
          col=c("lightblue","lightgreen"))
 
stripchart(list(pre,post),vertical=T,pch=16,method="jitter",cex=0.5,add=T)
 
for more see
 
surefoss page

Thursday, July 26, 2012

Thursday, June 14, 2012

Perl DBD::mysql on OS X Lion

I ran into a problem this week when the perl libraries for accessing MySQL databases refused to work on Mac OS X Lion.

The solution is

export DYLD_LIBRARY_PATH="/usr/local/mysql/lib:$DYLD_LIBRARY_PATH"

in your .bashrc 

Installing WGET for MAC OS X Lion (and others)


First, use curl to download the latest wget source:
curl -O http://ftp.gnu.org/gnu/wget/wget-1.13.4.tar.gz
Next we use tar to uncompress the files you just downloaded:
tar -xzf wget-1.13.4.tar.gz
Use cd to change to the directory:
cd wget-1.13.4
Configure with the appropriate –with-ssl flag to prevent a “GNUTLS not available” error:
./configure --with-ssl=openssl
Build the source:
make
Install wget, it ends up in /usr/local/bin/:
sudo make install
Confirm everything worked by running wget:
wget --help
Clean up by removing wget source files when finished:
cd .. && rm -rf wget*
You’re all set, enjoy wget in Mac OS X.

Monday, December 26, 2011

About Reference Human Genome Assembly

Modernizing Reference Genome Assemblies

Church DM, Schneider VA, Graves T, Auger K, Cunningham F, et al. (2011)
Modernizing Reference Genome Assemblies.
PLoS Biol 9(7): e1001091.
doi:10.1371/journal.pbio.1001091

A good way to understand current standard of Assembly curation.

Friday, March 11, 2011

Another R-package to draw heatmap

pheatmap package on CRAN allows you to draw nice heatmaps - especially it allows you to control the cell size.

pheatmap(matrix, scale="row", cluster_cols=FALSE, cluster_rows=FALSE, cellwidth=10, cellheight=10, fontsize=8)

Tuesday, September 07, 2010

Generating random numbers in R

Original link : http://blog.revolutionanalytics.com/2009/02/how-to-choose-a-random-number-in-r.html


Generate a random number between 5.0 and 7.5

If you want to generate a decimal number where any value (including fractional values) between the stated minimum and maximum is equally likely, use the runif function. This function generates values from the Uniform distribution. Here's how to generate one random number between 5.0 and 7.5:

> x1 <- runif(1, 5.0, 7.5)
> x1
[1] 6.715697

Of course, when you run this, you'll get a different number, but it will definitely be between 5.0 and 7.5. You won't get the values 5.0 or 7.5 exactly, either.

If you want to generate multiple random values, don't use a loop. You can generate several values at once by specifying the number of values you want as the first argument to runif. Here's how to generate 10 values between 5.0 and 7.5:

> x2 <- runif(10, 5.0, 7.5)
> x2
[1] 6.339188 5.311788 7.099009 5.746380 6.720383 7.433535 7.159988
[8] 5.047628 7.011670 7.030854

Generate a random integer between 1 and 10

This looks like the same exercise as the last one, but now we only want whole numbers, not fractional values. For that, we use the sample function:

> x3 <- sample(1:10, 1)
> x3
[1] 4

The first argument is a vector of valid numbers to generate (here, the numbers 1 to 10), and the second argument indicates one number should be returned. If we want to generate more than one random number, we have to add an additional argument to indicate that repeats are allowed:

> x4 <- sample(1:10, 5, replace=T)
> x4
[1] 6 9 7 6 5

Note the number 6 appears twice in the 5 numbers generated. (Here's a fun exercise: what is the probability of running this command and having no repeats in the 5 numbers generated?)

Select 6 random numbers between 1 and 40, without replacement

If you wanted to simulate the lotto game common to many countries, where you randomly select 6 balls from 40 (each labelled with a number from 1 to 40), you'd again use the sample function, but this time without replacement:

> x5 <- sample(1:40, 6, replace=F)
> x5
[1] 10 21 29 12 7 31

You'll get a different 6 numbers when you run this, but they'll all be between 1 and 40 (inclusive), and no number will repeat. Also, you don't actually need to include the replace=F option -- sampling without replacement is the default -- but it doesn't hurt to include it for clarity.

Select 10 items from a list of 50

You can use this same idea to generate a random subset of any vector, even one that doesn't contain numbers. For example, to select 10 distinct states of the US at random:

> sample(state.name, 10)
[1] "Virginia" "Oklahoma" "Maryland" "Michigan"
[5] "Alaska" "South Dakota" "Minnesota" "Idaho"
[9] "Indiana" "Connecticut"

You can't sample more values than you have without allowing replacements:

> sample(state.name, 52)
Error in sample(state.name, 52) :
cannot take a sample larger than the population when 'replace = FALSE'

... but sampling exactly the number you do have is a great way to randomize the order of a vector. Here are the 50 states of the US, in random order:

> sample(state.name, 50)
[1] "California" "Iowa" "Hawaii"
[4] "Montana" "South Dakota" "North Dakota"
[7] "Louisiana" "Maine" "Maryland"
[10] "New Hampshire" "Rhode Island" "Texas"
[13] "Florida" "North Carolina" "Minnesota"
[16] "Arkansas" "Pennsylvania" "Colorado"
[19] "Idaho" "Connecticut" "Utah"
[22] "South Carolina" "Illinois" "Ohio"
[25] "New Jersey" "Indiana" "Wisconsin"
[28] "Mississippi" "Michigan" "Wyoming"
[31] "West Virginia" "Alaska" "Georgia"
[34] "Vermont" "Virginia" "Oklahoma"
[37] "Washington" "New Mexico" "New York"
[40] "Delaware" "Nevada" "Alabama"
[43] "Kentucky" "Missouri" "Oregon"
[46] "Tennessee" "Arizona" "Massachusetts"
[49] "Kansas" "Nebraska"

You could also have just used sample(state.name) for the same result -- sampling as many values as provided is the default.

Further reading

For more information about how R generates random numbers, check out the following help pages:

> ?runif
> ?sample
> ?.Random.seed

The last of these provides technical detail on the random number generator R uses, and how you can set the random seed to recreate strings of random numbers.

Tuesday, July 20, 2010

Running Geneset enrichment analysis on commandline

java -cp /stor1/shah/Ruben_peptide/gsea_analysis/gsea2-2.06.jar
-Xmx2000m xtools.gsea.Gsea
-res foo_expression_values.gct
-cls foo_expression_values.cls#Normal_versus_Treatment
-gmx msigdb.v2.5.symbols.gmt -chip HG_U133_Plus_2.chip
-collapse true -mode Max_probe -norm meandiv -nperm 1000
-permute phenotype -rnd_type no_balance -scoring_scheme weighted
-rpt_label my_analysis -metric Signal2Noise -sort real -order descending -include_only_symbols true -make_sets true -median false -num 100
-plot_top_x 20 -rnd_seed timestamp -save_rnd_lists false -set_max 500 -set_min 15 -zip_report false -out /stor1/shah/Ruben_peptide/gsea_analysis -gui false

commandline for GSEA.

Monday, July 19, 2010

Significance of overlapping gene lists

Wen Fury and Wentian Li
http://www.nslij-genetics.org/wli/pub/ieee-embs06.pdf

To identify significance of overlap for two differentially expressed gene sets n1 and n2 (e.g. d1-n1 and d2-n1) use either hypergeometric or Fisher's exact test p-value.

Given integers n, n1, n2, m (max(n1,n2) <= n and m <= min (n1,n2)), the hypergeometric distribution is defined as

P(m) = [C(n1, m) * C (n - n1, n2 -m)]/ C (n, n2)

where C(n,m) is the number of possibilities of choosing m objects out of n objects : C (n,m) = n!/[m! (n -m)!]

It is usually more interesting to calculate the sum of P(m) for m's equal or larger than the observed value (i.e. p-value) :


p-value = Sigma [k= m to min (n1,n2)] p(k)
= Sigma [k = 0 to min (n1,n2)] p(k) - Sigma [k = 0 to m - 1] p(k)

For calculating it in R use :

if m = 0, p-value = 1

phyper (m, n1, n - n1, n2):
p-value = phyper(min(n1,n2), n1, n-n1, n2) - phyper(m-1, n1, n-n1, n2) if m > 0

One can also use Fisher's exact test on the following 2-by-2 table:

col1 col2 total
row1 m n1-m n1
row2 n2-m n-n1-n2+m n - n1
total n2 n-n2 n

They produce identical results.

Thursday, July 08, 2010

From a logical matrix to numerical matrix

Que : From a matrix of TRUE/FALSE get a matrix of 0 and 1
Ans : Multiply the logical matrix by self

Monday, April 26, 2010

Removing all NA rows and columns

Removing all NA rows and/or columns

fsFit[-which(apply(fsFit,1,function(x)all(is.na(x)))),-which(apply(fsFit,2,function(x)all(is.na(x))))]

Wednesday, April 21, 2010

extracting a percentage of data by random by groups

1) Randomly choose 10% of data from each "age" group.

> x <- data.frame(group=sample(1:4,100,TRUE), age=runif(100,4,80))
> tapply(x$age, x$group, function(z) mean(z[sample(seq_along(z), length(z) / 10)]))


2) To split my dataset randomly into 2 parts: a prediction set (with 2/3 of my data) and a validation set (with 1/3 of my data).

> x <- 1:100 # test data
> y <- split(x, sample(1:2, length(x), replace=TRUE, prob=c(1,2)))

3) I would like to randomly divide this data frame in half. how to select those rows that were not selected and assign them to randomsample2

selected<-rep(0,39622)
selected[sample(1:39622,39622/2)]<-1
data$selected<-selected
rm(selected)
or
data$selected<-rbinom(39622,1,.5)