Skip to content

Commit ad502b9

Browse files
keiranmrainekeiranmraine
committed
Merge branch 'hotfix/v6.1.1'
2 parents 5e5c3c2 + 63824e9 commit ad502b9

File tree

10 files changed

+127
-82
lines changed

10 files changed

+127
-82
lines changed

‎CHANGES.md‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,11 @@
11
# Changes
22

3+
## v6.1.1
4+
5+
* Reduce I/O for small cpu overhead in coverage step
6+
* Fix assembly step to cope with bam.csi and cram.crai
7+
* R lib issues and disable diagnostic plots.
8+
39
## v6.1.0
410

511
* Add travis-ci

‎Rsupport/libInstall.R‎

Lines changed: 17 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -9,19 +9,27 @@ source("http://bioconductor.org/biocLite.R")
99
ipak <- function(pkg){
1010
new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
1111
if (length(new.pkg))
12-
biocLite(new.pkg, ask=FALSE, lib=instLib)
12+
biocLite(new.pkg, ask=FALSE, lib=instLib, lib.loc=instLib)
1313
sapply(pkg, library, character.only = TRUE)
1414
}
1515

16-
biocPackages <- c("data.table", "gam")
17-
ipak(biocPackages)
18-
19-
install.packages("VGAM_1.0-3.tar.gz", type="source", lib=instLib)
20-
21-
biocPackages <- c("stringr", "poweRlaw", "zlibbioc", "RColorBrewer")
22-
ipak(biocPackages)
23-
2416
install.packages("devtools", lib=instLib)
2517
library(devtools)
2618
options(download.file.method = "auto")
19+
20+
ipak(c("data.table"))
21+
ipak(c("gam"))
22+
23+
if ( version$major == 3 && version$minor < 2 ) {
24+
install.packages("VGAM_1.0-3.tar.gz", type="source", lib=instLib, lib.loc=instLib)
25+
} else {
26+
ipak(c("VGAM"))
27+
}
28+
29+
ipak(c("stringr"))
30+
ipak(c("mgcv"))
31+
ipak(c("poweRlaw"))
32+
ipak(c("zlibbioc"))
33+
ipak(c("RColorBrewer"))
34+
2735
install_github("sb43/copynumber", ref="f1688edc154f1a0e3aacf7781090afe02882f623")

‎perl/bin/brass-assemble‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -82,7 +82,7 @@ sub show_usage {
8282
my ($option) = @_;
8383

8484
my $usage = <<"EOF";
85-
Usage: brass-assemble [OPTION]... BEDPE-FILE [BAM-FILE]...
85+
Usage: brass-assemble [OPTION]... BEDPE-FILE [BAM|CRAM-FILE]...
8686
Options:
8787
-o FILE Write events to FILE rather than standard output
8888
-O FORMAT Write the next output stream in the specified FORMAT

‎perl/bin/compute_coverage.pl‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -45,8 +45,8 @@
4545

4646
use PCAP::Bam;
4747

48-
const my $PROPER_INTERSECT => q{bash -c 'set -o pipefail; samtools view -ub -f 2 -F 3852 -q 1 %s %s | bedtools intersect -loj -a %s -b stdin -sorted -c > %s'};
49-
const my $FOLDBK_INTERSECT => q{bash -c 'set -o pipefail; (samtools view -H %s; samtools view -F 3854 %s %s | %s %s) | samtools view -Sbu - | bedtools intersect -loj -a %s -b stdin -sorted -c > %s'};
48+
const my $PROPER_INTERSECT => q{bash -c 'set -o pipefail; samtools view -ub -f 2 -F 3852 -q 1 %s %s | bedtools intersect -loj -a %s -b stdin -sorted -c | bgzip -c > %s'};
49+
const my $FOLDBK_INTERSECT => q{bash -c 'set -o pipefail; (samtools view -H %s; samtools view -F 3854 %s %s | %s %s) | samtools view -Sbu - | bedtools intersect -loj -a %s -b stdin -sorted -c | bgzip -c > %s'};
5050

5151
if(scalar @ARGV < 3) {
5252
die "USAGE: ./compute_coverage.pl gc_windows.bed[.gz] in.bam out_path [chr_name]\n";
@@ -72,8 +72,8 @@
7272
# now need to capture all windows from original window file by the chr_name
7373
$new_windows = "$out/windows_$chr_name.bed";
7474
run_ext(qq{zgrep -wE '^$chr_name' $windows > $new_windows});
75-
$cn_bed_file = qq{$out/${sample_name}.$chr_name.ngscn.bed};
76-
$cn_fb_file = qq{$out/${sample_name}.$chr_name.ngscn.fb_reads.bed};
75+
$cn_bed_file = qq{$out/${sample_name}.$chr_name.ngscn.bed.gz};
76+
$cn_fb_file = qq{$out/${sample_name}.$chr_name.ngscn.fb_reads.bed.gz};
7777
}
7878
run_ext(sprintf $PROPER_INTERSECT, $bam, $chr_name, ($new_windows || $windows), $cn_bed_file);
7979
run_ext(sprintf $FOLDBK_INTERSECT, $bam, $bam, $chr_name, $^X, _which('brass_foldback_reads.pl'), $cn_bed_file, $cn_fb_file);

‎perl/bin/coverage_merge.pl‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -39,8 +39,8 @@
3939
use Const::Fast qw(const);
4040
use List::Util qw(first);
4141

42-
const my $FILE_PAIR => '%s.%s.ngscn.bed';
43-
const my $FILE_FOLD => '%s.%s.ngscn.fb_reads.bed';
42+
const my $FILE_PAIR => '%s.%s.ngscn.bed.gz';
43+
const my $FILE_FOLD => '%s.%s.ngscn.fb_reads.bed.gz';
4444

4545
die "Usage: genome.fa.fai sample_name indir include_list" unless(scalar @ARGV >= 3);
4646

@@ -90,7 +90,7 @@ sub cat_to_gzip {
9090
push @args, sprintf $format, $sample, $chr;
9191
die "Expected file missing $indir/$args[-1]\n" unless(-e $args[-1]);
9292
}
93-
my $command = qq{bash -c 'set -o pipefail; cat @args | gzip -c > $outfile'};
93+
my $command = qq{bash -c 'set -o pipefail; zcat @args | gzip -c > $outfile'};
9494
warn $command."\n";
9595
system($command) == 0 or die "Failed to merge files to $outfile: $!\n";
9696
}

‎perl/docs.tar.gz‎

1000 Bytes
Binary file not shown.

‎perl/lib/Bio/Brass.pm‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -46,7 +46,7 @@ our @EXPORT = qw(find_breakpoints find_dusty_vertices
4646
is_dusty get_isolated_bp_alignment get_isolated_bp_surrounding_region
4747
$VERSION);
4848

49-
our $VERSION = '6.1.0';
49+
our $VERSION = '6.1.1';
5050

5151
=head1 NAME
5252

‎perl/lib/Bio/Brass/ReadSelection.pm‎

Lines changed: 50 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -68,8 +68,8 @@ use constant { PROPER_PAIRED => 0x2, UNMAPPED => 0x4, MATE_UNMAP => 0x8,
6868
6969
$bam = Bio::Brass::ReadSelection->new($filename);
7070
71-
Opens the specified BAM file and associated BAI index file.
72-
$filename is either the BAM filename or "BAM:BAI" to also give the BAI index
71+
Opens the specified BAM|CRAM file and associated BAI|CSI|CRAI index file.
72+
$filename is either the alignment filename or "ALIGNMENTS:INDEX" to also give the index
7373
filename, in which case temporary symlinks will be created as necessary.
7474
7575
=cut
@@ -90,15 +90,33 @@ sub new {
9090
# Index file is .bam.bai, so Bio::DB::HTS can use the files directly
9191
undef $bam; undef $bai;
9292
}
93+
elsif (-e $bam && -e "$bam.csi") {
94+
# Index file is .bam.csi, so Bio::DB::HTS can use the files directly
95+
undef $bam; undef $bai;
96+
}
97+
elsif (-e $bam && $bam =~ m/cram$/ && -e "$bam.crai") {
98+
# Alignment is cram and index file is .cram.crai, so Bio::DB::HTS can use the files directly
99+
undef $bam; undef $bai;
100+
}
93101
else {
94-
$bai = $bam; $bai =~ s{[.][^./]*$}{}; $bai .= '.bai';
95-
if (-e $bam && -e $bai) {
96-
# Index file is .bai, so Bio::DB::HTS will need symlinks
97-
}
98-
else {
99-
# No index file present, so just fail with the original filename
100-
undef $bam; undef $bai;
101-
}
102+
# handle files which don't include bam/cram before index extension
103+
$bai = $bam;
104+
$bai =~ s{[.][^./]*$}{};
105+
if($bam =~ m/bam$/) {
106+
$bai .= '.bai';
107+
# doesn't matter if changed and not existing due to logic below
108+
$bai .= '.csi' unless(-e $bai);
109+
}
110+
elsif($bam =~ m/cram$/) {
111+
$bai .= '.crai';
112+
}
113+
if (-e $bam && -e $bai) {
114+
# Index file is .bai|.csi|.crai, so Bio::DB::HTS will need symlinks
115+
}
116+
else {
117+
# No index file present, so just fail with the original filename
118+
undef $bam; undef $bai;
119+
}
102120
}
103121

104122
if (defined $bam) {
@@ -109,21 +127,31 @@ sub new {
109127
# This is an abuse of tmpnam(), which tries to ensure that "$base" does
110128
# not exist, but says nothing about "$base.bam" or "$base.bam.bai"
111129
my $base = tmpnam();
130+
my $aln_ext = 'bam';
131+
$aln_ext = 'cram' if($bam =~ m/cram$/);
132+
my $idx_ext = 'bam.bai';
133+
if($bai =~ m/csi$/) {
134+
$idx_ext = 'bam.csi';
135+
}
136+
elsif($bai =~ m/cram$/) {
137+
$idx_ext = 'cram.crai';
138+
}
139+
112140

113-
symlink $bam, "$base.bam"
114-
or croak "can't symlink $bam to $base.bam: $!";
115-
unless (symlink $bai, "$base.bam.bai") {
116-
unlink "$base.bam"; # Ignore errors during destruction
117-
croak "can't symlink $bai to $base.bam.bai: $!";
141+
symlink $bam, "$base.$aln_ext"
142+
or croak "can't symlink $bam to $base.$aln_ext: $!";
143+
unless (symlink $bai, "$base.$idx_ext") {
144+
unlink "$base.$aln_ext"; # Ignore errors during destruction
145+
croak "can't symlink $bai to $base.$idx_ext: $!";
118146
}
119-
unless (symlink $bas, "$base.bam.bas") {
120-
unlink "$base.bam"; # Ignore errors during destruction
121-
unlink "$base.bam.bai"; # Ignore errors during destruction
122-
croak "can't symlink $bas to $base.bam.bas: $!";
147+
unless (symlink $bas, "$base.$aln_ext.bas") {
148+
unlink "$base.$aln_ext"; # Ignore errors during destruction
149+
unlink "$base.$idx_ext"; # Ignore errors during destruction
150+
croak "can't symlink $bas to $base.$aln_ext.bas: $!";
123151
}
124152

125153
$self->{unlink_base} = $base;
126-
$filename = "$base.bam";
154+
$filename = "$base.$aln_ext";
127155
}
128156

129157
$self->{bam} = Bio::DB::HTS->new(-bam => $filename);
@@ -175,7 +203,8 @@ sub DESTROY {
175203
my ($self) = @_;
176204
if (exists $self->{unlink_base}) {
177205
my $base = $self->{unlink_base};
178-
unlink "$base.bam", "$base.bam.bai"; # Ignore errors during destruction
206+
# Ignore errors during destruction
207+
unlink "$base.bam", "$base.bam.bai", "$base.bam.csi", "$base.cram", "$base.cram.crai";
179208
}
180209
}
181210

‎perl/lib/Sanger/CGP/Brass/Implement.pm‎

Lines changed: 7 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -90,8 +90,8 @@ const my $BRASS_FILTER => q{ -seq_depth 25.1 -blat %s -ref %s -tumour %s -infile
9090
#perl ~kr2/git/brass/perl/bin/brassI_filter.pl
9191

9292
## assemble
93-
const my $BRASS_ASSEMBLE => q{ -X -m mem -O bedpe -r %s -T %s -o %s %s %s:%s.bai %s:%s.bai};
94-
# extreme depth, genome.fa, tmp, output.tab, groups, tumourbam, tumourbam, normalbam, normalbam
93+
const my $BRASS_ASSEMBLE => q{ -X -m mem -O bedpe -r %s -T %s -o %s %s %s %s};
94+
# extreme depth, genome.fa, tmp, output.tab, groups, tumourbam, normalbam
9595

9696
## grass
9797
const my $GRASS => q{ -genome_cache %s -ref %s -species %s -assembly %s -platform %s -protocol %s -tumour %s -normal %s -file %s -add_header 'brassVersion=%s'};
@@ -233,7 +233,9 @@ sub normcn {
233233

234234
move($normcn_stub.'.abs_cn.bg',$r_stub.'.ngscn.abs_cn.bg') || die "Move failed: $!\n";
235235
move($normcn_stub.'.segments.abs_cn.bg', $r_stub.'.ngscn.segments.abs_cn.bg') || die "Move failed: $!\n";
236-
move($normcn_stub.'.diagnostic_plots.pdf', $r_stub.'.ngscn.diagnostic_plots.pdf') || die "Move failed: $!\n";
236+
if(-e $normcn_stub.'.diagnostic_plots.pdf') {
237+
move($normcn_stub.'.diagnostic_plots.pdf', $r_stub.'.ngscn.diagnostic_plots.pdf') || die "Move failed: $!\n";
238+
}
237239

238240
PCAP::Threaded::touch_success(File::Spec->catdir($tmp, 'progress'), 0);
239241
}
@@ -590,8 +592,8 @@ sub assemble {
590592
$tmp_assemble,
591593
$assembled,
592594
$split_file,
593-
$options->{'tumour'}, $options->{'tumour'},
594-
$options->{'normal'}, $options->{'normal'};
595+
$options->{'tumour'},
596+
$options->{'normal'};
595597

596598
PCAP::Threaded::external_process_handler(File::Spec->catdir($tmp, 'logs'), $command, $index);
597599

‎perl/share/Rscripts/normalise_cn_by_gc_and_fb_reads.R‎

Lines changed: 38 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -86,44 +86,44 @@ while (i < nrow(out_table)) {
8686
}
8787

8888
# skipping diagnostinc plots for chromosome 9
89-
90-
if(is.element("9",chrs)){
91-
cat("Creating diagnostic plots...\n", file = stderr())
92-
# Some diagnostic plots
93-
pdf(paste0(output_body, ".diagnostic_plots.pdf"))
94-
smoothScatter(gc, logratio, pch = ".", xlab = "GC-content", ylab = "Log-ratio")
95-
smoothScatter(tumour_fb_ratio, logratio, pch = ".", xlab = "Tumour FB ratio", ylab = "Log-ratio")
96-
smoothScatter(normal_fb_ratio, logratio, pch = ".", xlab = "Normal FB ratio", ylab = "Log-ratio")
97-
plot.gam(gam_m, rugplot = F)
98-
idx = d.t[,1] == "9"
99-
plot(d.t[idx, 3], d.t[idx, 6], xlab = "Chr 9 position", ylab = "Read depth", pch = ".")
100-
plot(d.t[idx, 3], normalised_logratio[idx], xlab = "Chr 9 position", ylab = "Normalised log-ratio", pch = ".")
101-
dev.off()
102-
}else if(is.element("chr9",chrs)){
103-
cat("Creating diagnostic plots...\n", file = stderr())
104-
# Some diagnostic plots
105-
pdf(paste0(output_body, ".diagnostic_plots.pdf"))
106-
smoothScatter(gc, logratio, pch = ".", xlab = "GC-content", ylab = "Log-ratio")
107-
smoothScatter(tumour_fb_ratio, logratio, pch = ".", xlab = "Tumour FB ratio", ylab = "Log-ratio")
108-
smoothScatter(normal_fb_ratio, logratio, pch = ".", xlab = "Normal FB ratio", ylab = "Log-ratio")
109-
plot.gam(gam_m, rugplot = F)
110-
idx = d.t[,1] == "chr9"
111-
plot(d.t[idx, 3], d.t[idx, 6], xlab = "Chr 9 position", ylab = "Read depth", pch = ".")
112-
plot(d.t[idx, 3], normalised_logratio[idx], xlab = "Chr 9 position", ylab = "Normalised log-ratio", pch = ".")
113-
dev.off()
114-
}else{
115-
cat("Diagnostic plots created for fifth chromosome in the list...\n", file = stderr())
116-
# Some diagnostic plots
117-
pdf(paste0(output_body, ".diagnostic_plots.pdf"))
118-
smoothScatter(gc, logratio, pch = ".", xlab = "GC-content", ylab = "Log-ratio")
119-
smoothScatter(tumour_fb_ratio, logratio, pch = ".", xlab = "Tumour FB ratio", ylab = "Log-ratio")
120-
smoothScatter(normal_fb_ratio, logratio, pch = ".", xlab = "Normal FB ratio", ylab = "Log-ratio")
121-
plot.gam(gam_m, rugplot = F)
122-
idx = d.t[,1] == chrs[5]
123-
plot(d.t[idx, 3], d.t[idx, 6], xlab = paste0("chr ",chrs[5]," position"), ylab = "Read depth", pch = ".")
124-
plot(d.t[idx, 3], normalised_logratio[idx], xlab = paste0("chr ",chrs[5]," position"), ylab = "Normalised log-ratio", pch = ".")
125-
dev.off()
126-
}
89+
library(mgcv)
90+
# if(is.element("9",chrs)){
91+
# cat("Creating diagnostic plots...\n", file = stderr())
92+
# # Some diagnostic plots
93+
# pdf(paste0(output_body, ".diagnostic_plots.pdf"))
94+
# smoothScatter(gc, logratio, pch = ".", xlab = "GC-content", ylab = "Log-ratio")
95+
# smoothScatter(tumour_fb_ratio, logratio, pch = ".", xlab = "Tumour FB ratio", ylab = "Log-ratio")
96+
# smoothScatter(normal_fb_ratio, logratio, pch = ".", xlab = "Normal FB ratio", ylab = "Log-ratio")
97+
# plot.gam(gam_m, rugplot = F)
98+
# idx = d.t[,1] == "9"
99+
# plot(d.t[idx, 3], d.t[idx, 6], xlab = "Chr 9 position", ylab = "Read depth", pch = ".")
100+
# plot(d.t[idx, 3], normalised_logratio[idx], xlab = "Chr 9 position", ylab = "Normalised log-ratio", pch = ".")
101+
# dev.off()
102+
# }else if(is.element("chr9",chrs)){
103+
# cat("Creating diagnostic plots...\n", file = stderr())
104+
# # Some diagnostic plots
105+
# pdf(paste0(output_body, ".diagnostic_plots.pdf"))
106+
# smoothScatter(gc, logratio, pch = ".", xlab = "GC-content", ylab = "Log-ratio")
107+
# smoothScatter(tumour_fb_ratio, logratio, pch = ".", xlab = "Tumour FB ratio", ylab = "Log-ratio")
108+
# smoothScatter(normal_fb_ratio, logratio, pch = ".", xlab = "Normal FB ratio", ylab = "Log-ratio")
109+
# plot.gam(gam_m, rugplot = F)
110+
# idx = d.t[,1] == "chr9"
111+
# plot(d.t[idx, 3], d.t[idx, 6], xlab = "Chr 9 position", ylab = "Read depth", pch = ".")
112+
# plot(d.t[idx, 3], normalised_logratio[idx], xlab = "Chr 9 position", ylab = "Normalised log-ratio", pch = ".")
113+
# dev.off()
114+
# }else{
115+
# cat("Diagnostic plots created for fifth chromosome in the list...\n", file = stderr())
116+
# # Some diagnostic plots
117+
# pdf(paste0(output_body, ".diagnostic_plots.pdf"))
118+
# smoothScatter(gc, logratio, pch = ".", xlab = "GC-content", ylab = "Log-ratio")
119+
# smoothScatter(tumour_fb_ratio, logratio, pch = ".", xlab = "Tumour FB ratio", ylab = "Log-ratio")
120+
# smoothScatter(normal_fb_ratio, logratio, pch = ".", xlab = "Normal FB ratio", ylab = "Log-ratio")
121+
# plot.gam(gam_m, rugplot = F)
122+
# idx = d.t[,1] == chrs[5]
123+
# plot(d.t[idx, 3], d.t[idx, 6], xlab = paste0("chr ",chrs[5]," position"), ylab = "Read depth", pch = ".")
124+
# plot(d.t[idx, 3], normalised_logratio[idx], xlab = paste0("chr ",chrs[5]," position"), ylab = "Normalised log-ratio", pch = ".")
125+
# dev.off()
126+
# }
127127

128128
cat("Outputting data...\n", file = stderr())
129129

0 commit comments

Comments
 (0)