Skip to content

Instantly share code, notes, and snippets.

import pandas as pd
import numpy as np
from scipy.signal import find_peaks
from collections import OrderedDict
test = "GMM_logit_pvalue"
df = pd.read_csv("out_nanocompore_results.tsv", sep="\t")
df["Peak"] = 0
df = df[["pos", "chr", "genomicPos", "ref_id", "strand", "ref_kmer", "Peak", test]]
transcripts = set(df["ref_id"])
@tleonardi
tleonardi / combine_pvalues_hou.py
Last active October 17, 2019 14:39
Hou's method for the approximation for the distribution of the weighted combination of non-independent or independent probabilities.
from scipy.stats import chi2
import numpy as np
def combine_pvalues_hou(pvalues, weights, cor_mat):
""" Hou's method for the approximation for the distribution of the weighted
combination of non-independent or independent probabilities.
If any pvalue is nan, returns nan.
https://doi.org/10.1016/j.spl.2004.11.028
pvalues: list of pvalues to be combined
weights: the weights of the pvalues
@tleonardi
tleonardi / plotDispEsts2.R
Created August 29, 2018 13:15
Replacement for DESeq2 plotDispEsts() tha uses ggplot2
plotDispEsts2 <- function(dds){
require(dplyr)
require(reshape2)
require(ggplot2)
as.data.frame(mcols(dds)) %>%
select(baseMean, dispGeneEst, dispFit, dispersion) %>%
melt(id.vars="baseMean") %>%
filter(baseMean>0) %>%
ggplot(aes(x=baseMean, y=value, colour=variable)) +
geom_point(size=0.1) +
@tleonardi
tleonardi / david2html
Created January 27, 2017 11:36
Converts DAVID functional annotations to html
# Converts a tab-separated file containing
# DAVID's functional annotation results to
# Markdown format
# It reports: KEGG and GO terms along with
# the gene symbol and name
FILE=
cut -f 1,2,5,6,7,9 $FILE | awk '
BEGIN{OFS=FS="\t"}NR>1{
print "#"$1
#!/bin/awk -f
function logfactorial(n, f, i){
f=0
for(i=1;i<=n;i++){
f+=log(i)
}
return f
}
# m: Mean
@tleonardi
tleonardi / tund.sh
Last active April 12, 2016 10:12
Function to make an ssh tunnel that forwards a local display port back to the login node's display port in an LSF cluster that doesn't support X Forwarding
# Tunnel Display
# This function sets up an ssh tunnel that forwards a local display port
# back to the login node's display port. The tunnel is controlled by the
# socket in $HOME/tmp/ssh-${LSB_SUB_HOST}-${HOSTNAME}-${DISPLAY_NUMBER}.
# The ssh tunnel is in the background and keeps running even after the
# interactive shell is closed, thus preventing completion of the LSF job.
# To avoid this, we setup a trap on SIGINT SIGTERM EXIT that uses the ssh
# control socket to signal the tunnel to exit.
tund(){
if [[ -n $LSB_JOBID ]]; then
@tleonardi
tleonardi / tophatScoresAndFlags.md
Created February 25, 2016 11:40
Snippets from TopHat source code that print the flags and scores of the alignments

#Flags and scores in Tophat output

Tophat version 2.1.1

There are two functions 'bool rewrite_sam_record()', with different sets of parameters.

The first one is used for singletons the second for paired alignements:

Keybase proof

I hereby claim:

  • I am tleonardi on github.
  • I am tleonardi (https://keybase.io/tleonardi) on keybase.
  • I have a public key whose fingerprint is E6BA BA6A 0F8A 5B99 24B9 25E0 5FB4 E671 D50A 5A4E

To claim this, I am signing this object:

@tleonardi
tleonardi / knitMe.Rmd
Last active September 7, 2015 10:04
Knitr and pandoc-fignos
```{r set-options}
my_hook <- function(x, options) {
if (options$fig.show == 'animate') return(hook_plot_html(x, options))
"%n%" <- knitr:::"%n%"
base = opts_knit$get('base.url') %n% ''
cap = knitr:::.img.cap(options)
if (is.null(w <- options$out.width) & is.null(h <- options$out.height) &
is.null(s <- options$out.extra) & options$fig.align == 'default') {
return(sprintf('![%s](%s%s){#fig:%s} ', cap, base, knitr:::.upload.url(x), options$label))
@tleonardi
tleonardi / tanimoto.R
Created February 23, 2015 11:52
Tanimoto distance in R
tanimoto <- function(x, similarity=F) {
res<-sapply(x, function(x1){
sapply(x, function(x2) {i=length(which(x1 & x2)) / length(which(x1 | x2)); ifelse(is.na(i), 0, i)})
})
if(similarity==T) return(res)
else return(1-res)
}
x <- data.frame(Samp1=c(0,0,0,1,1,1,0,0,1), Samp2=c(1,1,1,1,1,1,0,0,1), Samp3=c(1,0,0,1,0,0,0,0,0), Samp4=c(1,1,1,1,1,1,1,1,1))
tanimoto(x)