Changes between Version 1 and Version 2 of SOP/scRNA-seq/tradeSeq


Ignore:
Timestamp:
10/27/20 15:27:49 (4 years ago)
Author:
ibarrasa
Comment:

--

Legend:

Unmodified
Added
Removed
Modified
  • SOP/scRNA-seq/tradeSeq

    v1 v2  
    55[https://github.com/statOmics/tradeSeq tradeSeq GitHub page]
    66
    7 [https://github.com/statOmics/tradeSeq--- Input from From slingShot]
     7[https://bioconductor.org/packages/devel/bioc/vignettes/tradeSeq/inst/doc/tradeSeq.html#load-data tradeSeq tradeSeq workflow using input from From slingShot]
    88
    9 [https://github.com/statOmics/tradeSeq--- Input from From Monocle]
     9[https://statomics.github.io/tradeSeq/articles/Monocle.html Workflow using input from Monocle]
    1010
    11 Example code
     11Example code importing only counts pseudotime and cells weights on each trajectory
     12
     13{{{
     14#! /usr/bin/Rscript
     15
     16library(tradeSeq)
     17library(monocle)
     18library(RColorBrewer)
     19library(SingleCellExperiment)
     20library(slingshot)
     21RNGversion("3.5.0")
     22palette(brewer.pal(8, "Dark2"))
     23
     24packageDescription("tradeSeq")$Version
     25#load data
     26
     27cellPSMatrix <- as.matrix(read.table("Subset.PsTime.txt", sep="\t"))
     28cellWeights <- as.matrix(read.table("Subset.CellsWeight.txt", sep="\t"))
     29counts <- as.matrix(read.table("RawCounts.txt", sep="\t"))
     30
     31#this step is optional
     32#select the subset of cells that are present on the pseudotime and cells weight file
     33countsFiltered <- subset.matrix(counts, select=rownames(cellWeights)) 
     34
     35#Evaluate the knots that we'll use when fitting the General Additive Models
     36icMat <- evaluateK(counts = countsFiltered,  k = 3:10,
     37                   nGenes = 200, verbose = T,  pseudotime = cellPSMatrix,
     38                   cellWeights = cellWeights)
     39
     40library(BiocParallel)
     41library(magrittr)
     42
     43# Register BiocParallel Serial Execution (no parallelization in that case)
     44BiocParallel::register(BiocParallel::SerialParam())
     45FitModels <- fitGAM(counts = countsFiltered,nknots = 8, pseudotime = cellPSMatrix, cellWeights = cellWeights )
     46
     47patternRes <- patternTest(FitModels)
     48oPat <- order(patternRes$waldStat, decreasing = TRUE)
     49head(rownames(patternRes)[oPat])
     50save.image("Data_afterPatternTest.RData")
     51
     52pdf( "plotSmoothers.pdf", w=11, h=8.5)
     53plotSmoothers(FitModels, countsFiltered , gene = rownames(patternRes)[oPat][1])
     54plotSmoothers(FitModels, countsFiltered , gene = rownames(patternRes)[oPat][2])
     55plotSmoothers(FitModels, countsFiltered , gene = rownames(patternRes)[oPat][3])
     56plotSmoothers(FitModels, countsFiltered , gene = rownames(patternRes)[oPat][4])
     57plotSmoothers(FitModels, countsFiltered , gene = rownames(patternRes)[oPat][5])
     58dev.off()
     59
     60write.table(patternRes, file="PatternTest8knots.txt", sep="\t", quote=F, row.names=T)
     61
     62
     63}}}
     64