Additional examples

Additional usage examples

Analysis at the gene level with exon counts

Note: Be sure to have constructed a metaseqR2 annotation database prior to continuing with the following examples!

An additional example with human exon data follows. If you have a multiple core system, be careful on how you are using the restrictCores option and generally how many cores you are using with scripts purely written in R. The analysis with exon read data can cause memory problems when a lot of BAM files are imported, so unless you have more than 64Gb of RAM available, consider setting setting the above parameter to something like 0.2. Generally, although the metaseqR2 gene model construction process from exon data is much faster and memory-efficient now, it requires some caution especially for large experiments. As example BAM files from a realistic dataset that can demonstrate the full availabilities of metaseqR2 do not fit within the Bioconductor package, you can find additional examples small (3 chromosomes) example dataset here. The reference genome for this example is hg19. If you have already buit a metaseqR2 annotation database, please provide it with the localDb argument, otherwise, an on-the-fly annotation will be downloaded and constructed for this example.

# A full example pipeline with exon counts

# Download sample data
lnk <- "https://drive.google.com/uc?export=download&id=1AuFAJSGO3wwubH_Xq949ZHd_hB68wT3P"
download.file(lnk,file.path(tempdir(),"test_data.zip"))
unzip(file.path(tempdir(),"test_data.zip"),exdir=tempdir(),junkpaths=TRUE)

# Construct targets file
targetsDf <- data.frame(
    samplename=c("CON1","CON2","DOX1","DOX2"),
    filename=dir(tempdir(),pattern=".bam$",full.names=TRUE),
    condition=c("CON","CON","DOX","DOX"),
    paired=rep("single",4),
    stranded=rep("forward",4)
)
targetsFile <- file.path(tempdir(),"targets.txt")
write.table(targetsDf,file=targetsFile,sep="\t",quote=FALSE,row.names=FALSE)

# Run the pipeline
metaseqr2(
    sampleList=targetsFile,
    contrast="DOX_vs_CON",
    org="hg19",
    countType="exon",
    normalization="deseq",
    statistics="deseq",
    pcut=0.05,
    qcPlots=c(
         "mds","biodetection","countsbio","saturation","rnacomp","pairwise",
         "boxplot","gcbias","lengthbias","meandiff","meanvar","correl",
         "deheatmap","volcano","mastat","biodist","filtered"
    ),
    figFormat=c("png","pdf"),
    exportWhat=c("annotation","p_value","adj_p_value","fold_change",
        "stats","counts"),
    exportScale=c("natural","log2","log10","vst"),
    exportValues=c("raw","normalized"),
    exportStats=c("mean","median","sd","mad","cv","rcv"),
    restrictCores=0.5,
    geneFilters=list(
        length=list(
            length=500
        ),
        avgReads=list(
            averagePerBp=100,
            quantile=0.25
        ),
        expression=list(
            median=TRUE,
            mean=FALSE
        ),
        biotype=getDefaults("biotypeFilter","hg19")
    ),
    exportWhere=file.path(tempdir(),"test5")
)

or in a more simplified version

# A full example pipeline with exon counts
metaseqr2(
    sampleList=targetsFile,
    contrast="DOX_vs_CON",
    org="hg19",
    countType="exon",
    normalization="deseq",
    statistics="deseq",
    preset="medium_normal",
    restrictCores=0.5,
    exportWhere=file.path(tempdir(),"test6")
)

Analysis at the 3'UTR level

Note: As with the previous example, be sure to have constructed a metaseqR2 annotation database prior to continuing.

metaseqR2 includes also facilities to analyze data from protocols that sequence only the 3'UTR of a gene/transcript for gene expression profile monitoring. Such protocols include for example the Quant-Seq kit by Lexogen. You can find additional a downsampled Quant-Seq dataset mapped to mm10 here.

# A full example pipeline with 3'UTR counts

# Download sample data
lnk <- "https://drive.google.com/uc?export=download&id=1SMtTf3j4O_4ycERoHd9Tafatheno2Tah"

download.file(lnk,file.path(tempdir(),"test_data.zip"))
unzip(file.path(tempdir(),"test_data.zip"),exdir=tempdir(),junkpaths=TRUE)

# Construct targets file
targetsDf <- data.frame(
    samplename=c("WT6","WT7","KO4","KO5"),
    filename=dir(tempdir(),pattern=".bam$",full.names=TRUE),
    condition=c("WT","WT","KO","KO"),
    paired=rep("single",4),
    stranded=rep("forward",4)
)
targetsFile <- file.path(tempdir(),"targets.txt")
write.table(targetsDf,file=targetsFile,sep="\t",quote=FALSE,row.names=FALSE)

# Run the pipeline
metaseqr2(
    sampleList=targetsFile,
    contrast="KO_vs_WT",
    org="mm10",
    countType="utr",
    normalization="deseq2",
    statistics="deseq2",
    pcut=0.05,
    qcPlots=c(
         "mds","biodetection","countsbio","saturation","rnacomp","pairwise",
         "boxplot","gcbias","lengthbias","meandiff","meanvar","correl",
         "deheatmap","volcano","mastat","biodist","filtered"
    ),
    figFormat="png",
    exportWhat=c("annotation","p_value","adj_p_value","fold_change",
        "stats","counts"),
    exportScale=c("natural","log2","log10","vst"),
    exportValues=c("raw","normalized"),
    exportStats=c("mean","median","sd","mad","cv","rcv"),
    restrictCores=0.5,
    geneFilters=list(
        length=list(
            length=500
        ),
        avgReads=list(
            averagePerBp=100,
            quantile=0.25
        ),
        expression=list(
            median=TRUE,
            mean=FALSE
        ),
        biotype=getDefaults("biotypeFilter","mm10")
    ),
    exportWhere=file.path(tempdir(),"test6")
)

or in a more simplified version

# A full example pipeline with 3'UTR counts
metaseqr2(
    sampleList=targetsFile,
    contrast="KO_vs_WT",
    org="mm10",
    countType="utr",
    normalization="deseq",
    statistics="deseq",
    preset="medium_normal",
    restrictCores=0.5,
    exportWhere=file.path(tempdir(),"test6")
)