1+ use crate :: covcalc:: { parse_regions, Region } ;
2+ use crate :: filehandler:: { is_bed_or_gtf, read_bedfile} ;
3+ use crate :: filtering:: Alignmentfilters ;
14use pyo3:: prelude:: * ;
25use pyo3:: types:: PyList ;
36use rayon:: prelude:: * ;
47use rayon:: ThreadPoolBuilder ;
58use rust_htslib:: bam:: { self , Header , IndexedReader , Read , Reader , Writer } ;
6- use tempfile:: { Builder , TempPath } ;
79use std:: fs:: File ;
810use std:: io:: Write ;
9- use crate :: covcalc:: { parse_regions, Region } ;
10- use crate :: filtering:: Alignmentfilters ;
11- use crate :: filehandler:: { is_bed_or_gtf, read_bedfile} ;
11+ use tempfile:: { Builder , TempPath } ;
1212
1313#[ pyfunction]
1414pub fn r_alignmentsieve (
15- bamifile : & str , // input bamfile
16- ofile : & str , // output file
17- nproc : usize , // threads
18- filter_metrics : & str , // filter metrics file.
15+ py : Python ,
16+ bamifile : & str , // input bamfile
17+ ofile : & str , // output file
18+ nproc : usize , // threads
19+ filter_metrics : & str , // filter metrics file.
1920 filtered_out_readsfile : & str , // filtered_out_reads bam/bedfile.
20- verbose : bool , // verbose
21- shift : Py < PyList > , // python list of the shift to perform.
22- _bed : bool , // output format in BEDPE.
23- filterrnastrand : & str , // "forward", "reverse" or "None".
24- minmappingquality : u8 , // minimum mapping quality.
25- samflaginclude : u16 , // sam flag include
26- samflagexclude : u16 , // sam flag exclude
27- blacklist : & str , // blacklist file name.
28- minfraglen : u32 , // minimum fragment length.
29- maxfraglen : u32 , // maximum fragment length.
21+ verbose : bool , // verbose
22+ shift : Py < PyList > , // python list of the shift to perform.
23+ _bed : bool , // output format in BEDPE.
24+ filterrnastrand : & str , // "forward", "reverse" or "None".
25+ minmappingquality : u8 , // minimum mapping quality.
26+ samflaginclude : u16 , // sam flag include
27+ samflagexclude : u16 , // sam flag exclude
28+ blacklist : & str , // blacklist file name.
29+ minfraglen : u32 , // minimum fragment length.
30+ maxfraglen : u32 , // maximum fragment length.
3031 _extend_reads : u32 ,
3132 _center_reads : bool ,
32-
3333) -> PyResult < ( ) > {
3434 // Input bam file
3535 let bam = Reader :: from_path ( bamifile) . unwrap ( ) ;
@@ -40,13 +40,10 @@ pub fn r_alignmentsieve(
4040 if filtered_out_readsfile != "None" {
4141 write_filters = true ;
4242 }
43- let mut readshift: Vec < i32 > = Vec :: new ( ) ;
44- Python :: with_gil ( |py| {
45- readshift = shift. extract ( py) . expect ( "Failed to extract shift." ) ;
46- } ) ;
43+ let readshift: Vec < i32 > = shift. extract ( py) . expect ( "Failed to extract shift" ) ;
4744 // shift is of length 0, 2, or 4.
4845
49- // Define regions
46+ // Define regions
5047 let ( regions, chromsizes) = parse_regions ( "None" , vec ! [ bamifile] ) ;
5148 // If there is a blacklist, read it.
5249 let mut backlistregions: Option < Vec < Region > > = None ;
@@ -56,10 +53,11 @@ pub fn r_alignmentsieve(
5653 match isbed. as_str ( ) {
5754 "gtf" => panic ! ( "Error: Please provide a bed file for the blacklist." ) ,
5855 "bed" => {
59- let ( bls, _) = read_bedfile ( & blacklist. to_string ( ) , false , chromsizes. keys ( ) . collect ( ) ) ;
56+ let ( bls, _) =
57+ read_bedfile ( & blacklist. to_string ( ) , false , chromsizes. keys ( ) . collect ( ) ) ;
6058 backlistregions = Some ( bls) ;
61- } ,
62- _ => panic ! ( "Error: Cannot determine filetype of blacklist file." )
59+ }
60+ _ => panic ! ( "Error: Cannot determine filetype of blacklist file." ) ,
6361 }
6462 }
6563
@@ -79,17 +77,29 @@ pub fn r_alignmentsieve(
7977 ) ;
8078 let pool = ThreadPoolBuilder :: new ( ) . num_threads ( 1 ) . build ( ) . unwrap ( ) ;
8179 let ( sieve, filtersieve, totalreads, filteredreads) = pool. install ( || {
82- regions. par_iter ( )
83- . map ( |i| sieve_bamregion ( bamifile, i, & filters, & readshift, write_filters, nproc, verbose) )
80+ regions
81+ . par_iter ( )
82+ . map ( |i| {
83+ sieve_bamregion (
84+ bamifile,
85+ i,
86+ & filters,
87+ & readshift,
88+ write_filters,
89+ nproc,
90+ verbose,
91+ )
92+ } )
8493 . reduce (
8594 || ( Vec :: new ( ) , Vec :: new ( ) , 0 , 0 ) ,
86- |( mut _sieve, mut _filtersieve, mut _total, mut _filter) , ( sieve, filtersieve, total, filter) | {
95+ |( mut _sieve, mut _filtersieve, mut _total, mut _filter) ,
96+ ( sieve, filtersieve, total, filter) | {
8797 _sieve. extend ( sieve) ;
8898 _filtersieve. extend ( filtersieve) ;
8999 _total += total;
90100 _filter += filter;
91101 ( _sieve, _filtersieve, _total, _filter)
92- }
102+ } ,
93103 )
94104 } ) ;
95105
@@ -107,7 +117,8 @@ pub fn r_alignmentsieve(
107117 }
108118 // write filtered reads if necessary
109119 if write_filters {
110- let mut ofilterbam = Writer :: from_path ( filtered_out_readsfile, & header, bam:: Format :: Bam ) . unwrap ( ) ;
120+ let mut ofilterbam =
121+ Writer :: from_path ( filtered_out_readsfile, & header, bam:: Format :: Bam ) . unwrap ( ) ;
111122 let _ = ofilterbam. set_threads ( nproc) ;
112123 for sb in filtersieve. into_iter ( ) {
113124 if let Some ( sb) = sb {
@@ -127,15 +138,33 @@ pub fn r_alignmentsieve(
127138 // write header
128139 writeln ! ( of, "#bamFilterReads --filterMetrics" ) . unwrap ( ) ;
129140 writeln ! ( of, "#File\t Reads\t Remaining Total\t Initial Reads" ) . unwrap ( ) ;
130- writeln ! ( of, "{}\t {}\t {}" , bamifile, totalreads-filteredreads, totalreads) . unwrap ( ) ;
141+ writeln ! (
142+ of,
143+ "{}\t {}\t {}" ,
144+ bamifile,
145+ totalreads - filteredreads,
146+ totalreads
147+ )
148+ . unwrap ( ) ;
131149 }
132150
133151 Ok ( ( ) )
134152}
135153
136-
137- fn sieve_bamregion ( ibam : & str , regstruct : & Region , alfilters : & Alignmentfilters , _shift : & Vec < i32 > , write_filters : bool , nproc : usize , verbose : bool ) -> ( Vec < Option < TempPath > > , Vec < Option < TempPath > > , u64 , u64 ) {
138- let region = ( regstruct. chrom . clone ( ) , regstruct. get_startu ( ) , regstruct. get_endu ( ) ) ;
154+ fn sieve_bamregion (
155+ ibam : & str ,
156+ regstruct : & Region ,
157+ alfilters : & Alignmentfilters ,
158+ _shift : & Vec < i32 > ,
159+ write_filters : bool ,
160+ nproc : usize ,
161+ verbose : bool ,
162+ ) -> ( Vec < Option < TempPath > > , Vec < Option < TempPath > > , u64 , u64 ) {
163+ let region = (
164+ regstruct. chrom . clone ( ) ,
165+ regstruct. get_startu ( ) ,
166+ regstruct. get_endu ( ) ,
167+ ) ;
139168 let mut total_reads: u64 = 0 ;
140169 let mut filtered_reads: u64 = 0 ;
141170 let mut bam = IndexedReader :: from_path ( ibam) . unwrap ( ) ;
@@ -201,7 +230,7 @@ fn sieve_bamregion(ibam: &str, regstruct: &Region, alfilters: &Alignmentfilters,
201230 }
202231 continue ;
203232 }
204-
233+
205234 // SAM flags
206235 if alfilters. samflaginclude != 0 && ( record. flags ( ) & alfilters. samflaginclude ) == 0 {
207236 filtered_reads += 1 ;
@@ -223,7 +252,9 @@ fn sieve_bamregion(ibam: &str, regstruct: &Region, alfilters: &Alignmentfilters,
223252 // fragment length
224253 if alfilters. minfraglen != 0 || alfilters. maxfraglen != 0 {
225254 if record. is_paired ( ) {
226- if record. insert_size ( ) . abs ( ) < alfilters. minfraglen as i64 || record. insert_size ( ) . abs ( ) > alfilters. maxfraglen as i64 {
255+ if record. insert_size ( ) . abs ( ) < alfilters. minfraglen as i64
256+ || record. insert_size ( ) . abs ( ) > alfilters. maxfraglen as i64
257+ {
227258 filtered_reads += 1 ;
228259 if let Some ( filterbamout) = & mut filterbamout {
229260 filterbamout. write ( & record) . unwrap ( ) ;
@@ -264,7 +295,7 @@ fn sieve_bamregion(ibam: &str, regstruct: &Region, alfilters: &Alignmentfilters,
264295 }
265296 continue ;
266297 }
267- } ,
298+ }
268299 ( "forward" , false ) => {
269300 if !( record. flags ( ) & 16 == 16 ) {
270301 filtered_reads += 1 ;
@@ -274,7 +305,7 @@ fn sieve_bamregion(ibam: &str, regstruct: &Region, alfilters: &Alignmentfilters,
274305 }
275306 continue ;
276307 }
277- } ,
308+ }
278309 ( "reverse" , true ) => {
279310 if !( ( record. flags ( ) & 144 == 144 ) || ( record. flags ( ) & 96 == 96 ) ) {
280311 filtered_reads += 1 ;
@@ -284,7 +315,7 @@ fn sieve_bamregion(ibam: &str, regstruct: &Region, alfilters: &Alignmentfilters,
284315 }
285316 continue ;
286317 }
287- } ,
318+ }
288319 ( "reverse" , false ) => {
289320 if !( record. flags ( ) & 16 == 0 ) {
290321 filtered_reads += 1 ;
@@ -294,27 +325,33 @@ fn sieve_bamregion(ibam: &str, regstruct: &Region, alfilters: &Alignmentfilters,
294325 }
295326 continue ;
296327 }
297- } ,
298- _ => { } ,
328+ }
329+ _ => { }
299330 }
300331 }
301332 sievebamout. write ( & record) . unwrap ( ) ;
302333 written = true ;
303334 }
304335
305336 match ( written, filterwritten) {
306- ( true , true ) => {
307- ( vec ! [ Some ( sievebam_path) ] , vec ! [ Some ( filterbam_path) ] , total_reads, filtered_reads)
308- } ,
309- ( true , false ) => {
310- ( vec ! [ Some ( sievebam_path) ] , vec ! [ None ] , total_reads, filtered_reads)
311- } ,
312- ( false , true ) => {
313- ( vec ! [ None ] , vec ! [ Some ( filterbam_path) ] , total_reads, filtered_reads)
314- } ,
315- ( false , false ) => {
316- ( vec ! [ None ] , vec ! [ None ] , total_reads, filtered_reads)
317- }
337+ ( true , true ) => (
338+ vec ! [ Some ( sievebam_path) ] ,
339+ vec ! [ Some ( filterbam_path) ] ,
340+ total_reads,
341+ filtered_reads,
342+ ) ,
343+ ( true , false ) => (
344+ vec ! [ Some ( sievebam_path) ] ,
345+ vec ! [ None ] ,
346+ total_reads,
347+ filtered_reads,
348+ ) ,
349+ ( false , true ) => (
350+ vec ! [ None ] ,
351+ vec ! [ Some ( filterbam_path) ] ,
352+ total_reads,
353+ filtered_reads,
354+ ) ,
355+ ( false , false ) => ( vec ! [ None ] , vec ! [ None ] , total_reads, filtered_reads) ,
318356 }
319-
320- }
357+ }
0 commit comments