r2ascension/skill_repo/bio-epitranscriptomics-modification-visualization

library(Guitar)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
 
# Load m6A peaks
peaks <- import('m6a_peaks.bed')
 
# Create metagene plot
# Shows distribution relative to transcript features
GuitarPlot(
    peaks,
    txdb = TxDb.Hsapiens.UCSC.hg38.knownGene,
    saveToPDFprefix = 'm6a_metagene'
)

# Create bigWig from IP/Input ratio
bamCompare -b1 IP.bam -b2 Input.bam \
    --scaleFactors 1:1 \
    --ratio log2 \
    -o IP_over_Input.bw
 
# Metagene around stop codons
computeMatrix scale-regions \
    -S IP_over_Input.bw \
    -R genes.bed \
    --regionBodyLength 2000 \
    -a 500 -b 500 \
    -o matrix.gz
 
plotProfile -m matrix.gz -o metagene.pdf

# Create normalized bigWig for genome browser
bamCoverage -b IP.bam \
    --normalizeUsing CPM \
    -o IP_normalized.bw
 
# Peak BED to bigBed
bedToBigBed m6a_peaks.bed chrom.sizes m6a_peaks.bb

library(ComplexHeatmap)
 
# m6A signal around peaks
Heatmap(
    signal_matrix,
    name = 'm6A signal',
    cluster_rows = TRUE,
    show_row_names = FALSE
)