Skip to content

Instantly share code, notes, and snippets.

View ATpoint's full-sized avatar

Alexander Bender ATpoint

View GitHub Profile
library(recount)
library(edgeR)
#/ Get the counts from GTEx via recount as a SummarizedExperiment
#/ => 1.3GB file
options(timeout=600)
download_study("SRP012682", type = "rse-gene")
load(file.path("SRP012682", "rse_gene.Rdata"))
#/ remove whitespaces in tissue names:
$ head renv.lock -n 40
{
"R": {
"Version": "4.0.3",
"Repositories": [
{
"Name": "CRAN",
"URL": "https://cran.rstudio.com"
}
]
@ATpoint
ATpoint / stdout_problem.nf
Last active June 8, 2021 19:41
Problem: .nextflow.out piles up the entire stdout channel, feels like it should not do that when output is explicitely stdout
#!/usr/bin/env nextflow
nextflow.enable.dsl=2
// Dummy example: Decompress a BAM in processA => stdout => read in processB as stdin and compress back to BAM.
process processA {
input:
path(bam)
$ make
mkdir -p objs
g++ -std=c++11 -Wall -O3 -fopenmp -march=native -c src/chromap.cc -o objs/chromap.o -lm -lz
src/chromap.cc:1096:0: warning: ignoring #pragma omp taskloop [-Wunknown-pragmas]
#pragma omp taskloop grainsize(grain_size) //num_tasks(num_threads_* 50)
^
src/chromap.cc:2001:0: warning: ignoring #pragma omp taskloop [-Wunknown-pragmas]
#pragma omp taskloop num_tasks(num_threads_* num_threads_)
^
src/chromap.cc: In instantiation of 'void chromap::Chromap<MappingRecord>::MapSingleEndReads() [with MappingRecord = chromap::PAFMapping]':
library(DESeq2)
#/ from https://raw.githubusercontent.com/shangguandong1996/picture_link/main/WFX_count_Rmatrix.txt
data <- read.delim("WFX_count_Rmatrix.txt", row.names = "Geneid")
get_res <- function(dds){
dds <- DESeq(dds)
res <- results(dds, name=resultsNames(dds)[2])
lfcShrink(dds, coef=resultsNames(dds)[2], res=res, type="ashr")
}
@ATpoint
ATpoint / guides_angles.R
Last active September 24, 2021 10:35
No more messing with h/vjust when rotating labels
library(ggplot2)
library(reshape2)
library(egg)
# thanks to https://stackoverflow.com/questions/1330989/rotating-and-spacing-axis-labels-in-ggplot2/60650595#60650595
dat<-reshape2::melt(data.frame(groupA=rnorm(20),
groupB=rnorm(20),
groupC=rnorm(20)))
@ATpoint
ATpoint / getData.R
Last active December 18, 2022 10:06
# install.packages("data.table", "dplyr")
library(data.table)
library(dplyr)
parseProperly <- function(x){
x <- x[,-c(ncol(x)-1,ncol(x))]
x
}
@ATpoint
ATpoint / colorblind.R
Created December 19, 2022 22:39
A colorblind-friendly palette including gray and black
cb <- c("#E69F00", "#0072B2", "#009E73",
"#CC79A7", "#56B4E9", "#D55E00",
"#661100", "#F0E442", "#332288",
"gray40", "black")
library(ggplot2)
ggplot(data.frame(x=c(rep(c(1,2,3), each=3), 4,4), y=1:length(l), color=LETTERS[1:length(l)]),
aes(x=x, y=y, color=color)) +
geom_point(size=10, show.legend=FALSE) +
@ATpoint
ATpoint / convertUCell.R
Last active April 22, 2023 18:33
Convert UCell output to label
#' Take output of UCell::ScoreSignatures_UCell() and return a vector of labels with highest score.
#' If the highest score per row is unique (only one value is the highest value) return that label.
#' If all scores are zero return "NA".
#' If more than one value is the highest value return "ambiguous".
#' library(UCell)
#' data(sample.matrix)
#' gene.sets <- list(Tcell_signature = c("CD2","CD3E","CD3D"),
#' Myeloid_signature = c("SPI1","FCER1G","CSF1R"))
#' scores <- ScoreSignatures_UCell(sample.matrix, features=gene.sets)
#' convertUCell(scores)
@ATpoint
ATpoint / getPercentExpressed.R
Last active April 28, 2023 12:18
Calculate percentage of cells per group expressing the genes in a numeric matrix
# Given a matrix-like object with genes being rows and samples being columns,
# and a group information for the columns, calculate the percentage of columns per group
# that have counts > threshold per gene. Efficiently done via base R, compatible with dgCMatrix.
# data <- data.frame(A=c(1,1,0,0), B=c(1,0,0,0), C=c(1,0,0,0))
# group <- as.character(c(1,1,2))
get_pexpr <- function(data, group, threshold=0, digits=2){
if(ncol(data)!=length(group)) stop("ncol(data) != length(group)")
if(!is.numeric(threshold) | threshold < 0) stop("threshold must be numeric and > 0")