adamewing

samtools idxstats my.bam | awk '{reads += ($3+$4)} END {print reads}'

adamewing

#!/usr/bin/env python

''' Script for correcting broken mapsplice BAMs using pysam
    Adam Ewing (ewingad@soe.ucsc.edu)
'''

import sys
import pysam
import os
import subprocess

adamewing

#!/bin/sh

if [ $# -ne 1 ]
then
    cat << EOF

usage: $0 <list of SRA location, name (whitespace seperated)>

For example, make a file called sra.txt and add this to it:

adamewing

#!/usr/bin/env python

import vcf
import sys
import traceback
import gzip

'''
validatevcf.py
Contact: Adam Ewing (ewingad@soe.ucsc.edu)

adamewing

#!/usr/bin/env python

import pysam
import sys
import numpy as np

from itertools import izip
from random import uniform

def get_depth(bam, chrom, pos):

adamewing

#!/usr/bin/env python

import pysam
import sys
from re import sub

def getRG(tags):
    ''' fetch RG, tags is from a pysam.AlignedRead.tags, returns RG name '''
    for tag, val in tags:
        if tag == 'RG':

adamewing

#!/usr/bin/env python

import vcf
import sys

'''
PyVCF will parse INDELs as SVs if they have SVTYPE in the INFO field. 
This provides a quick-fix.

contact: adam.ewing@mater.uq.edu.au

adamewing

#!/usr/bin/env python

'''
Join and re-annotate repeatmasker entities.

Expects input to be from UCSC public mysql server repeatmasker tables, using a command similar to the following:
echo "select genoName, genoStart, genoEnd, repStart, repEnd, strand, repName from rmsk where repClass = 'Other'" | mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -D hg19 -N | sort -k1,1 -k2,2n

requires: pysam (for Fastafile), biopython, exonerate (https://www.ebi.ac.uk/~guy/exonerate/)

adamewing

#!/usr/bin/env python

'''
Simple copy number profiling script for tumour / normal pairs.

Adam Ewing
adam.ewing@mater.uq.edu.au

'''

adamewing

#!/usr/bin/env python

from collections import defaultdict as dd
from bx.intervals.intersection import Intersecter, Interval

def interval_forest(bed_file):
    ''' build dictionary of interval trees '''
    forest = dd(Intersecter)

    with open(bed_file, 'r') as bed: