Skip to content

Instantly share code, notes, and snippets.

slava ilnytskyy slavailn

  • Lethbridge
Block or report user

Report or block slavailn

Hide content and notifications from this user.

Learn more about blocking users

Contact Support about this user’s behavior.

Learn more about reporting abuse

Report abuse
View GitHub Profile
@slavailn
slavailn / SPIAcalc.R
Created Aug 12, 2019
Function to calculate SPIA scores for DESeq2 results table
View SPIAcalc.R
resFile <- "DESeq2_result_name" # has to have entrez ids as annotations
compName <- "comparison" # Name of the comparison
calcSPIA <- function(resFile, compName)
{
res <- read.table(resFile, header = T, sep = "\t", fill=T, quote = "\"", stringsAsFactors = F)
# Select differentially expressed genes (adjusted p-value < 0.05)
# This is a named vector, where names are entrez ids, and values are log2 fold changes
ids <- res[res$padj < 0.05,]$entrezgene_id
@slavailn
slavailn / collect_HISAT2_mapping_stats.sh
Created Jun 8, 2019
Collect HISAT2 mapping statistics from multiple summary files located in the current directory and organized them into a table
View collect_HISAT2_mapping_stats.sh
#! /bin/bash
# 22498713 reads; of these:
# 22498713 (100.00%) were unpaired; of these:
# 1754404 (7.80%) aligned 0 times
# 17667104 (78.52%) aligned exactly 1 time
# 3077205 (13.68%) aligned >1 times
# 92.20% overall alignment rate
echo -e "sampleID\ttotal_reads\tunmapped\taligned_one_time\taligned_multiple\talignment_rate"
@slavailn
slavailn / edgeR_GLM_smallRNA.R
Created Feb 20, 2019
Edger GLM approach with 2 treatments over 3 time-points on small RNA tags
View edgeR_GLM_smallRNA.R
library(edgeR)
setwd(<dir>)
files <- list.files("/tag_counts/", full.names = T) # tab delimited file with sequenca tags and raw counts
files
d <- readDGE(files, columns = c(1,2))
counts <- d$counts
sampleFile <- <sampleFile>
sampleInfo <- read.table(sampleFile, sep = "\t", header = T, stringsAsFactors = F)
@slavailn
slavailn / get_bwa_mapping_stats.sh
Created Jan 17, 2019
Get BWA mapping stats (mapped, unmapped reads, number of reads with MAPQ20)
View get_bwa_mapping_stats.sh
#! /bin/bash
echo -e "Samplename\tTotal_reads\tMapped_reads\tUnmapped_reads\tMAPQ20_reads\tPercent_mapped\tPercent_MAPQ20_reads"
for file in ./*bam
do
filename=`basename $file`
samplename=${filename%.bam}
total_reads=`samtools view -c $file`
mapped_reads=`samtools view -c -F 4 $file`
@slavailn
slavailn / Process_ChIP-seq_peaks.R
Created Jan 14, 2019
Process ChIP-seq peaks, based on multiple Bioconductor packages, but centered on ChipSeeker and ClusterProfiler
View Process_ChIP-seq_peaks.R
library(GenomicRanges)
library(ChIPpeakAnno)
library(genomation)
library(biomaRt)
library(ChIPseeker)
library(org.Mm.eg.db)
library(ReactomePA)
library(clusterProfiler)
setwd("<project_dir>")
@slavailn
slavailn / ChipSeeker_vignette.R
Created Jan 10, 2019
Going through ChIPseeker vignette, will modify as necessary for the actual analysis
View ChipSeeker_vignette.R
library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(clusterProfiler)
# Gene data from UCSC hg19
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
# Peak file in bed format (these are MACS files)
files <- getSampleFiles()
files
@slavailn
slavailn / GenomicDataCommons_rnaseq.R
Created Dec 21, 2018
GenomicDataCommons vignette examples and how to download RNA-seq data from GDC
View GenomicDataCommons_rnaseq.R
library(GenomicDataCommons)
library(DESeq2)
library(ggplot2)
setwd("~/Projects/MYT1L_TCGA_project/")
###############################################################################################
# Vignette examples
# Check connectivity and functionality
GenomicDataCommons::status()
@slavailn
slavailn / plotDMR_ggbio.R
Created Dec 12, 2018
plot DMR using ggbio (ggplot type graphics, separate pane for every sample with a bar graph for methylation percent)
View plotDMR_ggbio.R
library(GenomicRanges)
library(ggbio)
library(ChIPpeakAnno)
gr <- GRanges(seqnames = "NC_035082.1", IRanges(55355800, 55355870), strand = "*")
# Get methylation percents from metilene input and convert to GRanges object
input <- read.table("metilene_Control_Treated.input", sep = "\t", header = T, stringsAsFactors = F)
names(input)[1] <- "chr"
names(input)[2] <- "start"
@slavailn
slavailn / ggbio_examples.R
Created Dec 11, 2018
Going through code examples in ggbio
View ggbio_examples.R
library(GenomicRanges)
library(ggbio)
library(IRanges)
set.seed(1)
N <- 1000
gr <- GRanges(seqnames =
sample(c("chr1", "chr2", "chr3"),
size = N, replace = TRUE),
IRanges(
@slavailn
slavailn / show_genomic_tracks.R
Created Dec 10, 2018
Display genomic tracks - ggbio, sushi
View 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,
You can’t perform that action at this time.