Skip to content

Instantly share code, notes, and snippets.

View mtw's full-sized avatar

Michael T. Wolfinger mtw

46 followers · 65 following
View GitHub Profile
@mtw
mtw / GFF_change_strand
Last active December 20, 2015 19:29
Change strands in GFF file
cat ORIG.gff | perl -ne 'chomp; if($_ =~ /^\#/){print $_."\n";next;} my @F = split("\t"); if ($F[6] eq "+"){$F[6] = "-"}else{$F[6] = "+"};my $part1 = join("\t",$F[0],$F[1],$F[2],$F[3],$F[4],$F[5],$F[6],$F[7],$F[8]);print $part1."\n"' > ! AS.gff
@mtw
mtw / bam2fastqc.sh
Last active December 23, 2015 05:19
bam2fastqc
#!/bin/bash
fastqcdir="./FastQC"
if [ -d "$fastqcdir" ];
then
echo "$fastqcdir available"
else
mkdir $fastqcdir
fi
@mtw
mtw / count_mapped_reads.sh
Last active December 23, 2015 21:29
Wrapper script for counting mapped reads in a bam file
#!/bin/bash
countdir="./count"
count_mapped_reads=`which count_mapped_reads_fast.pl`
bam_stat=`which bam_stat.py`
samtools=`which samtools`
samstats=`which sam-stats`
if [ -d "$countdir" ];
then
echo "$countdir available"
@mtw
mtw / map_em_segemehl_pe.sh
Created October 21, 2013 11:24
Wrapper script for paired-end segemehl
#!/bin/bash
segemehl="/scratch/mtw/src/segemehl-svn/segemehl/segemehl.x"
results="./alignments"
reffa="../../genomes/TAIR10/TAIR10.fa"
refidx="../../genomes/TAIR10/TAIR10.idx"
samples=2
rep=3
if [ -d "$results" ];
@mtw
mtw / bam2bw.sh
Created October 22, 2013 11:39
Convert BAM to BigWig coverage (useful for visualization in UCSC Browser)
#!/bin/bash
# generate coverage profiles
vis="./vis"
chromsizes="foo.chrom.sizes"
genomeCoverageBed=`which genomeCoverageBed`
bedGraphToBigWig=`which bedGraphToBigWig`
if ! [ -d "$vis" ];
then
mkdir -p $vis
@mtw
mtw / htseq_count_em.sh
Last active December 26, 2015 08:19
Wrapper for htseq-count
#!/bin/bash
samdir="./"
annotation="my.gtf"
htseqdir="./htseq-count"
htseqcount=`which htseq-count`
samtools=`which samtools`
feature="gene"
idattr="gene_name"
ss="reverse"
@mtw
mtw / extract_pp.sh
Created October 26, 2013 22:05
Extract properly-paired reads and their mates (ie flags 99/147/163/83) from paired-end BAM files
#!/bin/bash
samtools=`which samtools`
bn=$(basename $1 .bam)
echo processing $bn
$samtools view -h -b -f99 $1 > $bn.flag99.bam
$samtools view -h -b -f147 $1 > $bn.flag147.bam
$samtools view -h -b -f163 $1 > $bn.flag163.bam
$samtools view -h -b -f83 $1 > $bn.flag83.bam
@mtw
mtw / cutadapt_em_PE.sh
Last active December 28, 2015 06:28
Wrapper script for paired-end cutadapt
#!/bin/bash
cutadapt=`which cutadapt`
fastqc=`which fastqc`
origdir=".."
results="."
fastqcdir="${results}/FastQC"
adapter5="CTACACTCTTTCCCTACACGACGCTCTTCCGATCT"
adapter3="GATCGGAAGAGCACACGTCTGAACTCCAGTCAC"
inprefix="C"
@mtw
mtw / dexseq_count_em.sh
Last active December 28, 2015 08:09
Wrapper for DEXseq read counting
#!/bin/bash
samdir="./"
dexseq_flat_gff="my.DEXSeq.gff"
outdir="./dexseq-count"
dexseq_count="python /home/foo/bin/dexseq_count.py"
samtools=`which samtools`
if ! [ -d "$outdir" ];
then
mkdir -p $outdir
@mtw
mtw / bam2fastq_cutadapt.sh
Last active February 7, 2017 22:56
Raw NGS data is often shipped as unaligned BAM files. This shell script converts raw reads to FASTQ format and trims Illumina adapters with cutadapt. Quality Control is performed with FastQC.
#!/bin/bash
cutadapt=`which cutadapt`
fastqc=`which fastqc`
bam2fastq=`which bam2fastq`
gzip=`which gzip`
origdir="."
results="cutadapt"
fastqcdir="${results}/FastQC"
adapter5="CTACACTCTTTCCCTACACGACGCTCTTCCGATCT"