jwaageSnippets

# substitute cut -f 20 with column containing sequence. Sed bit removes last line
cat file_with_sequences.bed | cut -f 20 | awk '{print; if (FNR % 1 == 0 ) printf ">" NR "\n";}' | sed -n '$!p' > sequences.fasta

jwaageSnippets

library("rtracklayer")
session <- browserSession("UCSC")
genome(session)<-"mm9"
query <- ucscTableQuery(session, "refGene")
tableName(query) <- "refGene"
getTable(query) -> refseq

jwaageSnippets

plot + scale_fill_discrete(name="Experimental\nCondition",
                         breaks=c("ctrl", "trt1", "trt2"),
                         labels=c("Control", "Treatment 1", "Treatment 2"))

jwaageSnippets

read.delim("/Volumes/BP_NextGen/jwg/Projects/Porse/TET1/homologene.data", stringsAsFactors=F) -> homologene
split(homologene, f=homologene[,1]) -> hSplit

#Tax ID: Homo 9606
#Mus: 10090

getHumanAndMouseHomologues <- function(x)
{
	if (c(9606)%in%x[,2]&&c(10090)%in%x[,2])
	{

jwaageSnippets

GTF2TranscriptDB <- function(gtf.file, out.file = NULL, verbose = TRUE)
{
    require(rtracklayer)
    require(GenomicRanges)
    require(GenomicFeatures)

    min.info <- c("gene_id", "transcript_id", "exon_number")

    if (verbose) message("Importing ", gtf.file)
    gtf <- import.gff(gtf.file, asRangedData = FALSE)

jwaageSnippets

sort -b -k18 inputfile.txt > outputfile.txt

jwaageSnippets

library("rtracklayer")
session <- browserSession("UCSC")
genome(session)<-"mm9"
query <- ucscTableQuery(session, "refGene")
tableName(query) <- "refGene"
getTable(query) -> refseq
refseq[,c(2,3,4,5,6,7,8,13)] -> refseq
refseq$"width" <- refseq$"txEnd"-refseq$"txStart"

as.character(refseq[,1]) -> refseq[,1]