Skip to content
Snippets Groups Projects
Commit aa2351e8 authored by Brandi Cantarel's avatar Brandi Cantarel
Browse files

updating statanal

parent 7adcb10c
Branches
Tags
No related merge requests found
......@@ -42,5 +42,5 @@ module load subread/1.6.1 stringtie/1.3.2d-intel
featureCounts -s $stranded -M --fraction -J --ignoreDup -T $SLURM_CPUS_ON_NODE -p -g gene_name -a ${gtf} -o ${pair_id}.cts ${sbam}
mkdir ${pair_id}_stringtie
cd ${pair_id}_stringtie
stringtie ../${sbam} -p $SLURM_CPUS_ON_NODE -G ../${gtf} -B -e -o denovo.gtf -A ../${pair_id}.fpkm.txt
stringtie ../${sbam} -p $SLURM_CPUS_ON_NODE -G ${gtf} -B -e -o denovo.gtf -A ../${pair_id}.fpkm.txt
#!/usr/bin/perl -w
#merge_tables.pl
use Getopt::Long qw(:config no_ignore_case no_auto_abbrev);
use File::Basename;
my %opt = ();
my $results = GetOptions (\%opt,'outdir|o=s','help|h');
$opt{outdir} = './' unless ($opt{outdir});
my @files = @ARGV;
foreach $file (@files) {
chomp($file);
open IN, "<$file" or die $!;
$fname = basename($file);
my $sample = (split(/\./,$fname))[0];
$samps{$sample} = 1;
my $head = <IN>;
chomp($head);
my @colnames = split(/\t/,$head);
while (my $line = <IN>) {
chomp($line);
my @row = split(/\t/,$line);
my $gene = $row[0];
my $ct = $row[-1];
$cts{$gene}{$sample} = $ct;
}
close IN;
}
open OUT, ">$opt{outdir}\/countFeature.txt" or die $!;
my @samples = sort {$a cmp $b} keys %samps;
print OUT join("\t",'FeatureCtType',@samples),"\n";
foreach $gene (keys %cts) {
my @line;
foreach $sample (@samples) {
unless ($cts{$gene}{$sample}) {
$cts{$gene}{$sample} = 0;
}
push @line, $cts{$gene}{$sample};
}
print OUT join("\t",$gene,@line),"\n";
}
close OUT;
#!/bin/bash
#rnaseqalign.sh
usage() {
echo "-h Help documentation for hisat.sh"
echo "-d --DEA"
echo "Example: bash statanal.sh"
exit 1
}
OPTIND=1 # Reset OPTIND
while getopts :d:h opt
do
case $opt in
d) dea=$OPTARG;;
h) usage;;
esac
done
shift $(($OPTIND -1))
# Check for mandatory options
if [[ -z $SLURM_CPUS_ON_NODE ]]
then
SLURM_CPUS_ON_NODE=1
fi
source /etc/profile.d/modules.sh
perl $baseDir/scripts/concat_cts.pl -o ./ *.cts
perl $baseDir/scripts/concat_fpkm.pl -o ./ *.fpkm.txt
perl $baseDir/scripts/concat_ctsum.pl -o ./ *.cts.summary
cp design.txt design.shiny.txt
cp geneset.gmt geneset.shiny.gmt
if [[ $dea == 'skip' ]]
then
touch empty.png
touch bg.rda
else
module load R/3.2.1-intel
Rscript $baseDir/scripts/dea.R
Rscript $baseDir/scripts/build_ballgown.R *_stringtie
perl $baseDir/scripts/concat_edgeR.pl *.edgeR.txt
fi
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment