extern crate crast;
extern crate num_cpus;
extern crate bzip2;
#[macro_use]
pub mod utls;
use std::thread;
use crast::crast;
use utls::*;
use bzip2::bufread::BzDecoder;
const DFLT_GP_OPN_PNLTY: AlgnScr = -7.;
const DFLT_GP_EXTNSN_PNLTY: AlgnScr = -1.;
const DFLT_THRSHLD_4_SD_FRQ: usize = 3;
const DFLT_THRSHLD_4_UNGP_ALGN_X_DRP: AlgnScr = 5.;
const DFLT_THRSHLD_4_GP_ALGN_X_DRP: AlgnScr = 15.;
const DFLT_CNTXT_MTCH_PRB: Prb = 0.25;
const DFLT_BS_RT: AlgnScr = 0.65;
const DFLT_THRSHLD_4_SD_STRCT_E_VL: Prb = 6.5e2;
const DFLT_THRSHLD_4_UNGP_ALGN_SQ_E_VL: Prb = 5.0e-3;
const DFLT_THRSHLD_4_UNGP_ALGN_STRCT_E_VL: Prb = 9.0e2;
const DFLT_THRSHLD_4_GP_ALGN_SQ_E_VL: Prb = 5.0e-4;
const DFLT_THRSHLD_4_GP_ALGN_STRCT_E_VL: Prb = 9.0e2;
const DFLT_THRSHLD_4_FNL_GP_ALGN_SQ_E_VL: Prb = 5.0e-4;
const DFLT_THRSHLD_4_FNL_GP_ALGN_STRCT_E_VL: Prb = 6.5e3;
fn main() {
let args = env::args().collect::<Vec<Arg>>();
let prgrm = args[0].clone();
let mut opts = Options::new();
opts.reqopt("", "db", "Database dir. of ref. seqs", "STR");
opts.reqopt("i", "", "Target RNA seqs in single file of FASTA format", "STR");
opts.reqopt("o", "", "Output file of MAF format", "STR");
opts.optopt("", "ms", &format!("Max. span between paired bases (Default: {})", DFLT_MX_SPN), "UINT");
opts.optopt("", "ts", &format!("Threshold of frequency of adaptive seeds. (Default: {})", DFLT_THRSHLD_4_SD_FRQ), "UINT");
opts.optopt("", "x1", &format!("Threshold of drop of score when ungapped alignment (Default: {})", DFLT_THRSHLD_4_UNGP_ALGN_X_DRP), "FLOAT");
opts.optopt("", "x2", &format!("Threshold of drop of score when gapped alignment (Default: {})", DFLT_THRSHLD_4_GP_ALGN_X_DRP), "FLOAT");
opts.optopt("", "cp", &format!("Match probability of secondary structure context (Default: {})", DFLT_CNTXT_MTCH_PRB), "FLOAT");
opts.optopt("b", "", &format!("Contribution ratio of base to score (Default: {})", DFLT_BS_RT), "FLOAT");
opts.optopt("", "e1", &format!("Threshold of # E-value of seeds with better scores on secondary structure per target seq. (Default: {:e})", DFLT_THRSHLD_4_SD_STRCT_E_VL), "FLOAT");
opts.optopt("", "e2", &format!("Threshold of # E-value of ungapped alignments with better scores on seq. per target seq. (Default: {:e})", DFLT_THRSHLD_4_UNGP_ALGN_SQ_E_VL), "FLOAT");
opts.optopt("", "e3", &format!("Threshold of # E-value of ungapped alignments with better scores on secondary structure per target seq. (Default: {:e})", DFLT_THRSHLD_4_UNGP_ALGN_STRCT_E_VL), "FLOAT");
opts.optopt("", "e4", &format!("Threshold of # E-value of gapped alignments with better scores on seq. per target seq. (Default: {:e})", DFLT_THRSHLD_4_GP_ALGN_SQ_E_VL), "FLOAT");
opts.optopt("", "e5", &format!("Threshold of # E-value of gapped alignments with better scores on secondary structure per target seq. (Default: {:e})", DFLT_THRSHLD_4_GP_ALGN_STRCT_E_VL), "FLOAT");
opts.optopt("", "e6", &format!("Threshold of # E-value of final gapped alignments with better scores on seq. per target seq. (Default: {:e})", DFLT_THRSHLD_4_GP_ALGN_SQ_E_VL), "FLOAT");
opts.optopt("", "e7", &format!("Threshold of # E-value of final gapped alignments with better scores on secondary structure per target seq. (Default: {:e})", DFLT_THRSHLD_4_FNL_GP_ALGN_STRCT_E_VL), "FLOAT");
opts.optopt("", "go", &format!("Gap opening penalty (Default: {})", DFLT_GP_OPN_PNLTY), "FLOAT");
opts.optopt("", "ge", &format!("Gap extension penalty (Default: {})", DFLT_GP_EXTNSN_PNLTY), "FLOAT");
opts.optopt("t", "", "# of threads in multithreading (Default: system val.)", "UINT");
opts.optopt("p", "", "Output file of target RNA context seqs for shortcut from next time (If not specified, don't produce)", "STR");
opts.optopt("u", "", "Input file of target RNA context seqs for shortcut (If not specified, don't use)", "STR");
opts.optflag("", "sg", "Alignment is semi-global (Default: false, local)");
opts.optflag("a", "", "Approximate estimation of target RNA context seqs");
opts.optflag("h", "help", "Print help menu");
let mtchs = match opts.parse(&args[1 ..]) {
Ok(mtch) => {mtch}
Err(fl) => {prnt_usg(&prgrm, &opts); panic!(fl.to_string())}
};
if mtchs.opt_present("h") {
prnt_usg(&prgrm, &opts);
return;
}
let db = mtchs.opt_str("db").expect("Failed to retrieve database dir. from command args");
let db = Path::new(&db);
let id_fl = BufReader::new(File::open(db.join(ID_FL)).expect("Failed to open seq. file"));
let mut id_fl = BzDecoder::new(id_fl);
let thrd_4_rf_sq_ids = thread::spawn(move || {
let mut bfr = String::new();
let _ = id_fl.read_to_string(&mut bfr);
bfr.lines().map(|id| {String::from(id)}).collect::<Vec<Id>>()
});
let sq_fl = BufReader::new(File::open(db.join(SQ_FL)).expect("Failed to open seq. file"));
let mut sq_fl = BzDecoder::new(sq_fl);
let thrd_4_rf_sqs = thread::spawn(move || {
let mut bfr = String::new();
let _ = sq_fl.read_to_string(&mut bfr);
bfr.lines().map(|sq| {Vec::from(sq.as_bytes())}).collect::<Vec<BSq>>()
});
let sfx_ar_fl = BufReader::new(File::open(db.join(SFX_AR_FL)).expect("Failed to open suffix array file"));
let mut sfx_ar_fl = BzDecoder::new(sfx_ar_fl);
let thrd_4_sfx_ars = thread::spawn(move || {
let mut bfr = String::new();
let _ = sfx_ar_fl.read_to_string(&mut bfr);
bfr.lines().map(|ln| {ln.split_whitespace().map(|chr| {chr.parse().expect("Failed to parse suffix index")}).collect()}).collect::<Vec<SfxAr>>()
});
let ptrn_hsh_mp_fl = BufReader::new(File::open(db.join(PTRN_HSH_MP_FL)).expect("Failed to open pattern hash map file"));
let mut ptrn_hsh_mp_fl = BzDecoder::new(ptrn_hsh_mp_fl);
let thrd_4_ptrn_hsh_mps = thread::spawn(move || {
let mut bfr = String::new();
let _ = ptrn_hsh_mp_fl.read_to_string(&mut bfr);
bfr.lines().map(|ln| {
let mut ptrn_hsh_mp = HashMap::default();
for sbstr in ln.split('\t') {
let elms = sbstr.split(',').collect::<Vec<&str>>();
if elms.len() == 3 {
ptrn_hsh_mp.insert(elms[0].split_whitespace().map(|belm| belm.parse().expect("Failed to parse element of bio seq.")).collect::<BSq>(), (elms[1].parse().expect("Failed to parse index of suffix array"), elms[2].parse().expect("Failed to parse index of suffix array")));
}
}
ptrn_hsh_mp
}).collect::<Vec<PtrnHshMp>>()
});
let i = mtchs.opt_str("i").expect("Failed to retrieve input file from command args");
let thrd_4_trgt_sq_itm_sts = thread::spawn(move || {
let mut trgt_sq_itm_sts: HashMap<_, _, Hshr> = HashMap::default();
let fasta_rdr = Reader::from_file(Path::new(&i)).expect("Failed to read FASTA file");
for rc in fasta_rdr.records() {
let dt = rc.expect("Failed to read FASTA record");
let sq = dt.seq();
trgt_sq_itm_sts.insert(String::from(dt.id().expect("Failed to get FASTA record ID")), (sq.to_vec(), Vec::new()));
}
trgt_sq_itm_sts
});
let t = if mtchs.opt_present("t") {
mtchs.opt_str("t").expect("Failed to retrieve # of threads from command args").parse().expect("Failed to parse # of threads")
} else {
num_cpus::get() as ThrdNm
};
let ts = if mtchs.opt_present("ts") {
mtchs.opt_str("ts").expect("Failed to retrieve threshold of frequency of adaptive seeds from command args").parse().expect("Failed to parse threshold of frequency of adaptive seeds")
} else {
DFLT_THRSHLD_4_SD_FRQ
};
let x1 = if mtchs.opt_present("x1") {
mtchs.opt_str("x1").expect("Failed to retrieve threshold of drop of score when ungapped alignment from command args").parse().expect("Failed to parse threshold of drop of score when ungapped alignment")
} else {
DFLT_THRSHLD_4_UNGP_ALGN_X_DRP
};
let x2 = if mtchs.opt_present("x2") {
mtchs.opt_str("x2").expect("Failed to retrieve threshold of drop of score when gapped alignment from command args").parse().expect("Failed to parse threshold of drop of score when gapped alignment")
} else {
DFLT_THRSHLD_4_GP_ALGN_X_DRP
};
let cp = if mtchs.opt_present("cp") {
mtchs.opt_str("cp").expect("Failed to retrieve match probability of RNA secondary structure context from command args").parse().expect("Failed to parse match probability of RNA secondary structure context")
} else {
DFLT_CNTXT_MTCH_PRB
};
let b = if mtchs.opt_present("b") {
mtchs.opt_str("b").expect("Failed to retrieve contribution ratio of base to score from command args").parse().expect("Failed to parse contribution ratio of base to score")
} else {
DFLT_BS_RT
};
let e1 = if mtchs.opt_present("e1") {
mtchs.opt_str("e1").expect("Failed to retrieve threshold of seed E-value on secondary structure from command args").parse().expect("Failed to parse threshold of seed E-value on secondary structure")
} else {
DFLT_THRSHLD_4_SD_STRCT_E_VL
};
let e2 = if mtchs.opt_present("e2") {
mtchs.opt_str("e2").expect("Failed to retrieve threshold of ungapped alignment E-value on seq. from command args").parse().expect("Failed to parse threshold of ungapped alignment E-value on seq.")
} else {
DFLT_THRSHLD_4_UNGP_ALGN_SQ_E_VL
};
let e3 = if mtchs.opt_present("e3") {
mtchs.opt_str("e3").expect("Failed to retrieve threshold of ungapped alignment E-value on secondary structure from command args").parse().expect("Failed to parse threshold of ungapped alignment E-value on secondary structure")
} else {
DFLT_THRSHLD_4_UNGP_ALGN_STRCT_E_VL
};
let e4 = if mtchs.opt_present("e4") {
mtchs.opt_str("e4").expect("Failed to retrieve threshold of gapped alignment E-value on seq. from command args").parse().expect("Failed to parse threshold of gapped alignment E-value on seq.")
} else {
DFLT_THRSHLD_4_GP_ALGN_SQ_E_VL
};
let e5 = if mtchs.opt_present("e5") {
mtchs.opt_str("e5").expect("Failed to retrieve threshold of gapped alignment E-value on secondary structure from command args").parse().expect("Failed to parse threshold of gapped alignment E-value on secondary structure")
} else {
DFLT_THRSHLD_4_GP_ALGN_STRCT_E_VL
};
let e6 = if mtchs.opt_present("e6") {
mtchs.opt_str("e6").expect("Failed to retrieve threshold of final gapped alignment E-value on seq. from command args").parse().expect("Failed to parse threshold of final gapped alignment E-value on seq.")
} else {
DFLT_THRSHLD_4_FNL_GP_ALGN_SQ_E_VL
};
let e7 = if mtchs.opt_present("e7") {
mtchs.opt_str("e7").expect("Failed to retrieve threshold of final gapped alignment E-value on secondary structure from command args").parse().expect("Failed to parse threshold of final gapped alignment E-value on secondary structure")
} else {
DFLT_THRSHLD_4_FNL_GP_ALGN_STRCT_E_VL
};
let ms = if mtchs.opt_present("ms") {
mtchs.opt_str("ms").expect("Failed to retrieve max. span between paired bases from command args").parse().expect("Failed to parse max. span between paired bases")
} else {
DFLT_MX_SPN
};
let go = if mtchs.opt_present("go") {
mtchs.opt_str("go").expect("Failed to retrieve gap opening penalty from command args").parse().expect("Failed to parse gap opening penalty")
} else {
DFLT_GP_OPN_PNLTY
};
let ge = if mtchs.opt_present("ge") {
mtchs.opt_str("ge").expect("Failed to retrieve gap opening penalty from command args").parse().expect("Failed to parse gap opening penalty")
} else {
DFLT_GP_EXTNSN_PNLTY
};
let prdcs_cntxt_dst_sqs;
let p = if mtchs.opt_present("p") {
prdcs_cntxt_dst_sqs = true;
mtchs.opt_str("p").expect("Failed to retrieve output file of target RNA context seqs from command args").parse().expect("Failed to parse output file of target RNA context seqs")
} else {
prdcs_cntxt_dst_sqs = false;
String::new()
};
let uss_cntxt_dst_sqs;
let u = if mtchs.opt_present("u") {
uss_cntxt_dst_sqs = true;
mtchs.opt_str("u").expect("Failed to retrieve input file of target RNA context seqs from command args").parse().expect("Failed to parse input file of target RNA context seqs")
} else {
uss_cntxt_dst_sqs = false;
String::new()
};
let o = mtchs.opt_str("o").expect("Failed to retrieve output file from command args");
let sg = mtchs.opt_present("sg");
let a = mtchs.opt_present("a");
let algn_scr_schm = (go, ge);
let mut bfr = String::new();
let mut cnfg_fl = BufReader::new(File::open(db.join(CNFG_FL)).expect("Failed to open config. file"));
let _ = cnfg_fl.read_to_string(&mut bfr);
let srchs_rvrs_strnd = !bfr.parse::<bool>().expect("Failed to parse if it aligns reverse seq.");
let algn_prms = (ts, x1, x2, cp, b, e1, (e2, e3), (e4, e5), (e6, e7));
let mut bfr = String::new();
let mut db_sz_fl = BufReader::new(File::open(db.join(DB_SZ_FL)).expect("Failed to open DB size file"));
let _ = db_sz_fl.read_to_string(&mut bfr);
let db_sz = bfr.parse::<usize>().expect("Failed to parse DB size");
let rf_sq_ids = thrd_4_rf_sq_ids.join().expect("Failed to get ref. seq. ids");
let rf_sq_nm = rf_sq_ids.len();
let prb_dst_sq_fl = BufReader::new(File::open(db.join(PRB_DST_SQ_FL)).expect("Failed to open prob. dist. seq. file"));
let mut prb_dst_sq_fl = BzDecoder::new(prb_dst_sq_fl);
let mut rf_sq_prb_dst_sqs = vec![Vec::new(); rf_sq_nm];
let mut thrd_pl = Pool::new(t);
let mut bfr = String::new();
let _ = prb_dst_sq_fl.read_to_string(&mut bfr);
thrd_pl.scoped(|scp| {
for (rf_sq_prb_dst_sq, lns) in rf_sq_prb_dst_sqs.iter_mut().zip(bfr.trim().split("\n\n")) {
scp.execute(move || {
*rf_sq_prb_dst_sq = lns.lines().map(|ln| {let prb_dst = ln.split_whitespace().map(|prb| prb.parse().expect("Failed to parse context prob.")).collect::<Vec<Prb>>(); [prb_dst[0], prb_dst[1], prb_dst[2], prb_dst[3], prb_dst[4], prb_dst[5]]}).collect();
});
}
});
let tmp_dr = String::from("/tmp/crast_") + >_tm_stmp();
let tmp_dr = Path::new(&tmp_dr);
if !tmp_dr.exists() {
let _ = create_dir(tmp_dr);
}
let tmp_cntxt_dst_sq_dr = tmp_dr.join("cntxt_dst_sqs");
let tmp_cntxt_dst_sq_dr = Path::new(&tmp_cntxt_dst_sq_dr);
if !tmp_cntxt_dst_sq_dr.exists() {
let _ = create_dir(tmp_cntxt_dst_sq_dr);
}
let mut trgt_sq_itm_sts = thrd_4_trgt_sq_itm_sts.join().expect("Failed to get target seq. item set"); // ID, seq., and context dist. seq.
if uss_cntxt_dst_sqs {
let mut cntxt_dst_sq_fl = BufReader::new(File::open(Path::new(&u)).expect("Failed to open input file"));
let mut bfr = String::new();
let _ = cntxt_dst_sq_fl.read_to_string(&mut bfr);
for lns in bfr.trim().split("\n\n") {
let mut lns = lns.lines();
let frst_ln = lns.next().expect("Failed to get first line");
trgt_sq_itm_sts.get_mut(&frst_ln[1 ..]).expect("Failed to get value of hash map").1 = lns.map(|ln| {let cntxt_dst = ln.split_whitespace().map(|prb| prb.parse().expect("Failed to parse context prob.")).collect::<Vec<Prb>>(); [cntxt_dst[0], cntxt_dst[1], cntxt_dst[2], cntxt_dst[3], cntxt_dst[4], cntxt_dst[5]]}).collect();
}
} else {
thrd_pl.scoped(|scp| {
for (ky, trgt_sq_itms) in trgt_sq_itm_sts.iter_mut() {
scp.execute(move || {
let cntxt_dst_sq = cpr(&trgt_sq_itms.0, ms, a);
trgt_sq_itms.1 = cntxt_dst_sq.iter().map(|cntxt_dst| [cntxt_dst[0] as Prb, cntxt_dst[1] as Prb, cntxt_dst[2] as Prb, cntxt_dst[3] as Prb, cntxt_dst[4] as Prb, cntxt_dst[5] as Prb]).collect();
if prdcs_cntxt_dst_sqs {
wrt_cntxt_dst_sq(ky, &cntxt_dst_sq, tmp_cntxt_dst_sq_dr);
}
});
}
});
}
if prdcs_cntxt_dst_sqs {
let mut otpt_fl_4_cntxt_dst_sqs = BufWriter::new(File::create(Path::new(&p)).expect("Failed to create output file"));
for cntxt_dst_sq in tmp_cntxt_dst_sq_dr.read_dir().expect("Failed to read entries within temp. dir.") {
let cntxt_dst_sq = cntxt_dst_sq.expect("Failed to read entry within temp. dir.");
let cntxt_dst_sq_pth = cntxt_dst_sq.path();
let mut cntxt_dst_sq = BufReader::new(File::open(cntxt_dst_sq_pth).expect("Failed to open entry within temp. dir."));
let mut bfr = Vec::new();
let _ = cntxt_dst_sq.read_to_end(&mut bfr);
let _ = otpt_fl_4_cntxt_dst_sqs.write_all(&bfr);
}
}
let tmp_pr_algn_dr = tmp_dr.join("pr_algns");
let tmp_pr_algn_dr = Path::new(&tmp_pr_algn_dr);
if !tmp_pr_algn_dr.exists() {
let _ = create_dir(tmp_pr_algn_dr);
}
let rf_sqs = thrd_4_rf_sqs.join().expect("Failed to get ref. seqs");
let sfx_ars = thrd_4_sfx_ars.join().expect("Failed to get suffix arrays");
let ptrn_hsh_mps = thrd_4_ptrn_hsh_mps.join().expect("Failed to get pattern hash maps");
if srchs_rvrs_strnd {
let rvrs_sfx_ar_fl = BufReader::new(File::open(db.join(RVRS_SFX_AR_FL)).expect("Failed to open suffix array file"));
let mut rvrs_sfx_ar_fl = BzDecoder::new(rvrs_sfx_ar_fl);
let mut bfr = String::new();
let _ = rvrs_sfx_ar_fl.read_to_string(&mut bfr);
let rvrs_sfx_ars = bfr.lines().map(|ln| {ln.split_whitespace().map(|chr| {chr.parse().expect("Failed to parse suffix index")}).collect()}).collect::<Vec<SfxAr>>();
let rvrs_ptrn_hsh_mp_fl = BufReader::new(File::open(db.join(RVRS_PTRN_HSH_MP_FL)).expect("Failed to open pattern hash map file"));
let mut rvrs_ptrn_hsh_mp_fl = BzDecoder::new(rvrs_ptrn_hsh_mp_fl);
let mut bfr = String::new();
let _ = rvrs_ptrn_hsh_mp_fl.read_to_string(&mut bfr);
let rvrs_ptrn_hsh_mps = bfr.lines().map(|ln| {
let mut rvrs_ptrn_hsh_mp = HashMap::default();
for sbstr in ln.split('\t') {
let elms = sbstr.split(',').collect::<Vec<&str>>();
if elms.len() == 3 {
rvrs_ptrn_hsh_mp.insert(elms[0].split_whitespace().map(|belm| belm.parse().expect("Failed to parse element of bio seq.")).collect::<BSq>(), (elms[1].parse().expect("Failed to parse index of suffix array"), elms[2].parse().expect("Failed to parse index of suffix array")));
}
}
rvrs_ptrn_hsh_mp
}).collect::<Vec<PtrnHshMp>>();
let rvrs_prb_dst_sq_fl = BufReader::new(File::open(db.join(RVRS_PRB_DST_SQ_FL)).expect("Failed to open reverse prob. dist. seq. file"));
let mut rvrs_prb_dst_sq_fl = BzDecoder::new(rvrs_prb_dst_sq_fl);
let mut bfr = String::new();
let _ = rvrs_prb_dst_sq_fl.read_to_string(&mut bfr);
let rf_sq_rvrs_prb_dst_sqs = bfr.trim().split("\n\n").map(|lns| {lns.lines().map(|ln| {let prb_dst = ln.split_whitespace().map(|chr| {chr.parse().expect("Failed to parse context prob.")}).collect::<Vec<Prb>>(); [prb_dst[0], prb_dst[1], prb_dst[2], prb_dst[3], prb_dst[4], prb_dst[5]]}).collect()}).collect::<Vec<PrbDstSq>>();
thrd_pl.scoped(|scp| {
for (rf_sq_id, rf_sq, sfx_ar, ptrn_hsh_mp, rf_sq_prb_dst_sq, rvrs_sfx_ar, rvrs_ptrn_hsh_mp, rf_sq_rvrs_prb_dst_sq) in multizip((rf_sq_ids.iter(), rf_sqs.iter(), sfx_ars.iter(), ptrn_hsh_mps.iter(), rf_sq_prb_dst_sqs.iter(), rvrs_sfx_ars.iter(), rvrs_ptrn_hsh_mps.iter(), rf_sq_rvrs_prb_dst_sqs.iter())) {
for (ky, trgt_sq_itms) in &trgt_sq_itm_sts {
scp.execute(move || {
crast(&(rf_sq_id, ky), &(rf_sq, &trgt_sq_itms.0), &vec![sfx_ar, rvrs_sfx_ar], &vec![ptrn_hsh_mp, rvrs_ptrn_hsh_mp], &vec![rf_sq_prb_dst_sq, rf_sq_rvrs_prb_dst_sq, &trgt_sq_itms.1], &algn_scr_schm, &algn_prms, srchs_rvrs_strnd, sg, db_sz, rf_sq_nm, &tmp_pr_algn_dr);
});
}
}
});
} else {
thrd_pl.scoped(|scp| {
for (rf_sq_id, rf_sq, sfx_ar, ptrn_hsh_mp, rf_sq_prb_dst_sq) in multizip((rf_sq_ids.iter(), rf_sqs.iter(), sfx_ars.iter(), ptrn_hsh_mps.iter(), rf_sq_prb_dst_sqs.iter())) {
for (ky, trgt_sq_itms) in &trgt_sq_itm_sts {
scp.execute(move || {
crast(&(rf_sq_id, ky), &(rf_sq, &trgt_sq_itms.0), &vec![sfx_ar], &vec![ptrn_hsh_mp], &vec![rf_sq_prb_dst_sq, &trgt_sq_itms.1], &algn_scr_schm, &algn_prms, srchs_rvrs_strnd, sg, db_sz, rf_sq_nm, &tmp_pr_algn_dr);
});
}
}
});
}
let mut otpt_fl = BufWriter::new(File::create(Path::new(&o)).expect("Failed to create output file"));
let _ = otpt_fl.write_all(b"##maf\n# CRAST version 1.0.4\n");
let bfr = format!("# Gap opening penalty = {}, gap extension penalty = {}, threshold of frequency of adaptive seeds = {}, threshold of drop of score when ungapped alignment = {}, threshold of drop of score when gapped alignment = {}, match probability of secondary structure context = {}, contribution ratio of base to score = {}, threshold of # E-value of seeds with better scores on secondary structure per target seq. = {:e}, threshold of # E-value of ungapped alignments with better scores on seq. per target seq. = {:e}, threshold of # E-value of ungapped alignments with better scores on secondary structure per target seq. = {:e}, threshold of # E-value of gapped alignments with better scores on seq. per target seq. = {:e}, threshold of # E-value of gapped alignments with better scores on secondary structure per target seq. = {:e}, threshold of # E-value of final gapped alignments with better scores on seq. per target seq. = {:e}, threshold of # E-value of final gapped alignments with better scores on secondary structure per target seq. = {:e}, max. span between paired bases = {}, {} alignment\n", go, ge, ts, x1, x2, cp, b, e1, e2, e3, e4, e5, e6, e7, ms, if sg {"semi-global"} else {"local"});
let _ = otpt_fl.write_all(bfr.as_bytes());
let bfr = format!("# Ref. seqs = {}, database size = {}\n", rf_sq_nm, db_sz);
let _ = otpt_fl.write_all(bfr.as_bytes());
let bfr = format!("# Reported # E-values are of alignments with better scores across database\n#E1 = E-value for primary sequence identity and E2 = E-value for secondary structure identity\n");
let _ = otpt_fl.write_all(bfr.as_bytes());
let bfr = String::from("# Coordinates are 0-based\n#For - strand matches, coordinates in reverse complements of ref. seqs are used\n");
let _ = otpt_fl.write_all(bfr.as_bytes());
let bfr = String::from("# Format = {id} {start} {alignment size} {strand} {seq. size} {alignment}\n\n");
let _ = otpt_fl.write_all(bfr.as_bytes());
for pr_algns in tmp_pr_algn_dr.read_dir().expect("Failed to read entries within temp. dir.") {
let pr_algns = pr_algns.expect("Failed to read entry within temp. dir.");
let mut pr_algns = BufReader::new(File::open(pr_algns.path()).expect("Failed to open entry within temp. dir."));
let mut bfr = Vec::new();
let _ = pr_algns.read_to_end(&mut bfr);
let _ = otpt_fl.write_all(&bfr);
}
let _ = remove_dir_all(tmp_dr);
}