Display genomic tracks - ggbio, sushi

show_genomic_tracks.R

library(Rsamtools)
library(rtracklayer)
library(Sushi)

# Draw raw alignment track for each bam file
setwd("~/Projects/Wiseman_Trout_RRBS/")
list.files("bam")

# Create TxDb object from GFF
chromInfo <- read.table("~/GENOMES/Rainbow_trout/GCF_002163495.1_Omyk_1.0_genomic.fa.fai", header = F,
                        stringsAsFactors = F, sep = "\t")
head(chromInfo)
chromInfo <- chromInfo[,c(1,2)]
names(chromInfo) <- c("chrom", "length")
is_circular <- rep(FALSE, dim(chromInfo)[1])
chromInfo$is_circular <- is_circular

gff_file <- "~/GENOMES/Rainbow_trout/GCF_002163495.1_Omyk_1.0_genomic.gff"

metadata <- data.frame(name="Resource URL", 
                       value=paste0("ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/002/163/495/GCF_002163495.1_Omyk_1.0/GCF_002163495.1_Omyk_1.0_genomic.gff.gz"))

TxDb.OMykiss.Refseq.Omyk.10.knownGene <- makeTxDbFromGFF(file=gff_file,
                                                         chrominfo=chromInfo,
                                                         dataSource="NCBI Refseq",
                                                         organism="Oncorhynchus mykiss",
                                                         metadata=metadata)
TxDb.OMykiss.Refseq.Omyk.10.knownGene


gr <- GRanges("NC_035082.1", IRanges(55355800-50000, 55355870+50000))

JM10 <- readGAlignments("bam/10JM_trimmed_bismark_bt2.sorted.bam", 
                      use.names=TRUE, param=ScanBamParam(which=gr))

JM1 <- readGAlignments("bam/1JM_trimmed_bismark_bt2.sorted.bam", 
                        use.names=TRUE, param=ScanBamParam(which=gr))

JM2 <- readGAlignments("bam/2JM_trimmed_bismark_bt2.sorted.bam", 
                       use.names=TRUE, param=ScanBamParam(which=gr))

JM3 <- readGAlignments("bam/3JM_trimmed_bismark_bt2.sorted.bam", 
                       use.names=TRUE, param=ScanBamParam(which=gr))

JM4 <- readGAlignments("bam/4JM_trimmed_bismark_bt2.sorted.bam", 
                       use.names=TRUE, param=ScanBamParam(which=gr))

JM5 <- readGAlignments("bam/5JM_trimmed_bismark_bt2.sorted.bam", 
                       use.names=TRUE, param=ScanBamParam(which=gr))

JM6 <- readGAlignments("bam/6JM_trimmed_bismark_bt2.sorted.bam", 
                       use.names=TRUE, param=ScanBamParam(which=gr))

JM7 <- readGAlignments("bam/7JM_trimmed_bismark_bt2.sorted.bam", 
                       use.names=TRUE, param=ScanBamParam(which=gr))

JM8 <- readGAlignments("bam/8JM_trimmed_bismark_bt2.sorted.bam", 
                       use.names=TRUE, param=ScanBamParam(which=gr))

JM9 <- readGAlignments("bam/9JM_trimmed_bismark_bt2.sorted.bam", 
                       use.names=TRUE, param=ScanBamParam(which=gr))

# Plot raw alignments
p10 <- autoplot(JM10, geom = "bar")
p1 <- autoplot(JM1, geom = "bar")
p2 <- autoplot(JM2, geom = "bar")
p3 <- autoplot(JM3, geom = "bar")
p4 <- autoplot(JM4, geom = "bar")
p5 <- autoplot(JM5, geom = "bar")
p6 <- autoplot(JM6, geom = "bar")
p7 <- autoplot(JM7, geom = "bar")
p8 <- autoplot(JM8, geom = "bar")
p9 <- autoplot(JM9, geom = "bar")

# gene track
p11 <- autoplot(TxDb.OMykiss.Refseq.Omyk.10.knownGene, which=gr, names.expr = "gene_id")

tracks(JM1=p1, JM2=p2, JM3=p3, JM4=p4, JM5=p5, JM6=p6,
       JM7=p7, JM8=p8, JM9=p9, JM10=p10, transcripts=p11, heights = c(rep(0.3, 11))) + ylab("")

# Plot coverage
p10 <- autoplot(JM10, geom = "bar", stat = "coverage", fill = "darkgreen")
p1 <- autoplot(JM1, geom = "bar", stat = "coverage", fill = "darkgreen")
p2 <- autoplot(JM2, geom = "bar", stat = "coverage", fill = "red")
p3 <- autoplot(JM3, geom = "bar", stat = "coverage", fill = "red")
p4 <- autoplot(JM4, geom = "bar", stat = "coverage", fill = "darkgreen")
p5 <- autoplot(JM5, geom = "bar", stat = "coverage", fill = "darkgreen")
p6 <- autoplot(JM6, geom = "bar", stat = "coverage", fill = "darkgreen")
p7 <- autoplot(JM7, geom = "bar", stat = "coverage", fill = "red")
p8 <- autoplot(JM8, geom = "bar", stat = "coverage", fill = "darkgreen")
p9 <- autoplot(JM9, geom = "bar", stat = "coverage", fill = "darkgreen")

tracks(JM1=p1, JM2=p2, JM3=p3, JM4=p4, JM5=p5, JM6=p6,
       JM7=p7, JM8=p8, JM9=p9, JM10=p10, transcripts=p11, heights = c(rep(0.3, 11))) + ylab("")+ ylim(c(0, 20))

# Import filtered bedgraphs
JM10 <- read.table("filtered_bedgraph/10JM_trimmed_bismark_bt2.sorted.bedgraph", sep = "\t")
JM9 <- read.table("filtered_bedgraph/9JM_trimmed_bismark_bt2.sorted.bedgraph", sep = "\t")
JM8 <- read.table("filtered_bedgraph/8JM_trimmed_bismark_bt2.sorted.bedgraph", sep = "\t")
JM7 <- read.table("filtered_bedgraph/7JM_trimmed_bismark_bt2.sorted.bedgraph", sep = "\t")
JM6 <- read.table("filtered_bedgraph/6JM_trimmed_bismark_bt2.sorted.bedgraph", sep = "\t")
JM5 <- read.table("filtered_bedgraph/5JM_trimmed_bismark_bt2.sorted.bedgraph", sep = "\t")
JM4 <- read.table("filtered_bedgraph/4JM_trimmed_bismark_bt2.sorted.bedgraph", sep = "\t")
JM3 <- read.table("filtered_bedgraph/3JM_trimmed_bismark_bt2.sorted.bedgraph", sep = "\t")
JM2 <- read.table("filtered_bedgraph/2JM_trimmed_bismark_bt2.sorted.bedgraph", sep = "\t")
JM1 <- read.table("filtered_bedgraph/1JM_trimmed_bismark_bt2.sorted.bedgraph", sep = "\t")

tiff("filtered_coverage_track.tiff", width = 500, height = 900)
par(mfrow=c(10,1))
par(mar=c(1,1,1,1))
plotBedgraph(signal=JM1, chrom="NC_035082.1", chromstart=55355800-100000, chromend=55355870+100000, lwd = 1, main="1JM")
plotBedgraph(signal=JM2, chrom="NC_035082.1", chromstart=55355800-100000, chromend=55355870+100000, lwd = 1,
             color="red", main="2JM")
plotBedgraph(signal=JM3, chrom="NC_035082.1", chromstart=55355800-100000, chromend=55355870+100000, lwd = 1,
             color="red", main="3JM")
plotBedgraph(signal=JM4, chrom="NC_035082.1", chromstart=55355800-100000, chromend=55355870+100000, lwd = 1, main="4JM")
plotBedgraph(signal=JM5, chrom="NC_035082.1", chromstart=55355800-100000, chromend=55355870+100000, lwd = 1, main="5JM")
plotBedgraph(signal=JM6, chrom="NC_035082.1", chromstart=55355800-100000, chromend=55355870+100000, lwd = 1, main="6JM")
plotBedgraph(signal=JM7, chrom="NC_035082.1", chromstart=55355800-100000, chromend=55355870+100000, lwd = 1,
             color="red", main="7JM")
plotBedgraph(signal=JM8, chrom="NC_035082.1", chromstart=55355800-100000, chromend=55355870+100000, lwd = 1, main="8JM")
plotBedgraph(signal=JM9, chrom="NC_035082.1", chromstart=55355800-100000, chromend=55355870+100000, lwd = 1, main="9JM")
plotBedgraph(signal=JM10, chrom="NC_035082.1", chromstart=55355800-100000, chromend=55355870+100000, lwd = 1, main="10JM")
dev.off()