diff --git a/.gitignore b/.gitignore index ad48c69..cc724e8 100644 --- a/.gitignore +++ b/.gitignore @@ -106,6 +106,22 @@ ENV/ .idea/ tests/ +!tests/ +tests/* +!tests/fixtures/ +!tests/fixtures/** +!tests/test_mini_fixtures_manifest.py +!tests/test_split_merge_python_expected.py +!tests/test_parse_se_python.py +!tests/test_cwl_yaml_adapter.py +tests/__pycache__/ +tests/*.pyc + +.conda-env/ +.snakemake/ +results/ +tests/fixtures/mini/source/failed_jobs_list.txt +tests/fixtures/mini/source/failed_jobs_list.txt.lck .ipynb_checkpoints/ SOURCEME diff --git a/NOTES b/NOTES deleted file mode 100644 index 16ceb68..0000000 --- a/NOTES +++ /dev/null @@ -1,7 +0,0 @@ - - -perl versions: 5.10, 5.12, 5.14, 5.16 have deterministic hashes - -5.18 started haivng non-deterministic hashes - -TSCC system-wide perl is 5.10 diff --git a/README.md b/README.md index 41ca0ee..bc38fa8 100644 --- a/README.md +++ b/README.md @@ -1,111 +1,128 @@ # repetitive-element-mapping -pipeline for mapping repetitive elements - -# Requirements: -- Bowtie2=2.2.6 -- Tested with perl=5.22.0 and perl=5.10.1 (perl 5.18+ introduces some randomness in -iterating over hashes, which may yield slightly different results). This is -due to the way this pipeline accesses reads via hashes, and determines ties -in quality by selecting the first read it encounters at the de-duplication -step. We have tried to mitigate this randomness by iterating over sorted -hash keys in 0.0.2. See NOTES. -- cwlref-runner=1.0 -- python 2.7 with pandas and numpy installed - -# Installation: -- For Yeo Lab: ```module load ecliprepmap``` -- For all others: - - see the ```DockerRequirement``` from each tool to see the specific requirements for each step. The custom scripts have been tested using the following software/versions: - - python=3.6 - - perl (5.10.1+) - - bowtie2=2.2.6 - - pandas (1.1.3+) - - numpy (1.18+) - - cwltool==1.0.20180306140409 - - cwltest - - galaxy-lib==17.9.3 - - ensure that the ```bin/```, ```bin/perl/```, and ```wf/```, and ```cwl/``` - are correctly in your path. - -# Example data: -- [example reference data for GRCh38 (hg38)](https://external-collaborator-data.s3-us-west-1.amazonaws.com/reference-data/repeat-family-mapping-grch38.tar.gz) -- [example input files from eCLIP for GRCh38 (hg38)](https://external-collaborator-data.s3-us-west-1.amazonaws.com/reference-data/example_data_for_repeat_mapping_hg38.tar.gz) -# Methods: -- (```map_repetitive_elements_/*.cwl``` - ```parse_bowtie2_output_realtime_includemultifamily_/*.pl```) Runs bowtie2 using the following commands: + +Python + Snakemake pipeline for repetitive element mapping with deterministic mini fixtures for validation. + +## What is in this branch +- Python implementation scripts in `bin/python/`. +- Snakemake workflow entrypoint in `Snakefile`. +- Rule modules in `workflow/rules/`. +- Validation scripts in `workflow/scripts/`. +- Mini fixtures and expected outputs in `tests/fixtures/mini/`. + +## Repository layout +- `Snakefile`: top-level workflow entrypoint. +- `config/config.yaml`: runtime file paths for mini validation inputs and expected outputs. +- `bin/python/split_bam_to_subfiles_SEorPE.py`: split SAM/BAM by UMI prefix. +- `bin/python/merge_multiple_parsed_files.simplified_20191022.py`: merge parsed statistics files. +- `workflow/rules/se_foundation.smk`: currently implemented, testable workflow stages. +- `workflow/scripts/verify_split_manifest.py`: checksum validation for split outputs. +- `workflow/scripts/verify_merge_against_expected.py`: checksum validation for merged output. + +## Requirements +- Python 3.11+ +- Snakemake 9+ +- samtools 1.23+ +- bowtie2 2.5+ +- pytest 9+ + +## Installation + +### 1) Clone +```bash +git clone https://github.com/YeoLab/repetitive-element-mapping.git +cd repetitive-element-mapping +git checkout codex/python-conversion +``` + +### 2) Create environment +```bash +mamba create -y -p ./.conda-env -c conda-forge -c bioconda \ + python=3.11 snakemake bowtie2 samtools pandas numpy pyyaml pytest +``` + +### 3) Activate +```bash +mamba activate ./.conda-env +``` + +## Input specification + +## CWL-style YAML Input Alignment +This branch supports two equivalent ways to provide run configuration: + +1) Populate top-level keys in `config/config.yaml`. +2) Pass an existing CWL-style job YAML at runtime: +```bash +HOME=$(pwd) ./.conda-env/bin/snakemake -j1 -p all \ + --config cwl_input_yaml=/absolute/path/to/job.yaml +``` + +Supported CWL-compatible keys: +- `barcode1r1FastqGz`, `barcode1rmRepBam` +- `barcode1Inputr1FastqGz`, `barcode1InputrmRepBam` +- `bowtie2_db`, `bowtie2_prefix`, `fileListFile1`, `gencodeGTF`, `gencodeTableBrowser`, `repMaskBEDFile` +- `dataset`, `prefixes`, `se_or_pe` + +For `class: File` / `class: Directory` objects, the `path` value is extracted. Relative `path` values are resolved relative to the YAML file location. + +The separate `mini_validation` block is only for deterministic fixture tests while conversion is in progress. +See `docs/config_mapping.md` for the one-to-one mapping details. +The current implemented Snakemake stages use mini fixture inputs configured in `config/config.yaml`: + +- `mini_validation.source_sam_gz`: compressed SAM-like file used for split stage. +- `mini_validation.merge_input_1`: parsed stats input file #1 for merge stage. +- `mini_validation.merge_input_2`: parsed stats input file #2 for merge stage. +- `mini_validation.split_manifest`: expected checksums for 25 split output files. +- `mini_validation.merge_expected`: expected merged parsed output. + +You can update these paths to point to your own test data as long as formats match. + +## Output specification +Running the current workflow stages produces: + +- `results/mini/split/python_tmp/*.tmp`: 25 prefix-split files (`AA`..`NN`). +- `results/mini/split/verified.ok`: split-stage validation success marker. +- `results/mini/merge/merged.python.parsed`: merged parsed output. +- `results/mini/merge/verified.ok`: merge-stage validation success marker. + +## Run instructions (deployment) + +### Local execution +```bash +HOME=$(pwd) ./.conda-env/bin/snakemake -j1 -p all +``` + +`HOME=$(pwd)` keeps Snakemake cache files inside workspace and avoids host cache permission issues. + +### Re-run all stages from scratch +```bash +HOME=$(pwd) ./.conda-env/bin/snakemake -j1 -p -F all +``` + +### Cluster deployment pattern +Use your site profile/launcher as usual, for example: +```bash +HOME=$(pwd) ./.conda-env/bin/snakemake --profile all +``` + +## Testing + +### Unit/integration tests ```bash -bowtie2 -q --sensitive -a -p 3 --no-mixed --reorder -x $bowtie_db -1 $fastq_file1 -2 $fastq_file2 2> $bowtie_out +./.conda-env/bin/python -m pytest -q tests/test_mini_fixtures_manifest.py ``` -where: - - - $fastq_file1 is read 1 of a trimmed CLIPSEQ expt - - $fastq_file2 is read 2 of a trimmed CLIPSEQ expt - - $bowtie_db is a bowtie database created using a manually curated set of repeat elements from RepBase or other. - - For every proper read pair: - - Mismatch score equals the sum of mismatch scores (as determined by the AS: field) for both reads - - Reads are compared on the basis of their mismatch scores and read quality, keeping the best or prioritize if equal. - - priority is based on input repeat database (e.g. primary transcripts are prioritized over pseudogenes). It is based on the ordering of the repeat database file (MASTER_filelist.wrepbaseandtRNA.enst2id.fixed.UpdatedSimpleRepeat) - - Since one read may map to multiple elements, we must compare each element - - If read maps to multiple elements of the same family with equal score: - - Read mapping information (fastq line) is kept for the element with the highest priority. - - Name of all elements is kept - - Note: all reads that map to multiple families (not used for downstream analysis). - - Create a SAM-like file -- (```deduplicate.cwl``` - ```duplicate_removal_inline_paired.count_region_other_reads_masksnRNAs_andreparse_SEandPE_20201210_simple.pl```) Merge repeat analysis with unique genomic mapping and remove PCR duplicates - - Use the randomer UMIs to remove duplicates based on quality - - To save memory, both the SAM-like file and the uniquely mapped BAM file - are split based on the first two nucleotides of the randomner umi. De-duplication - is performed serially for AA, AC, AG, AT .. NN. - - Between a read mapping to both unique genomic and repeat family, if - unique mapping to genome is more than 2 mismatches per - read = 2 * 2 * 6 alignment score better than to repeat element, - throw out repeat element and use genome mapping. Otherwise keep the repeat - element mapped read. - -### Determining the most correct assignment among elements within one family: -- Between longer and shorter transcripts: keep the shorter one -- If a read maps to two places on the same transcript, keep the first -- Treat "rRNA extra hash" different. See: ```parse_bowtie2_output_realtime_includemultifamily.pl``` $rRNA_extra_hash - - This is different because the rRNA precusor transcript contains 18S, 28S, 5.8S transcripts which are treated as separate families. -# Outputs: -.parsed file: tabbed file containing 4 comment lines and 4 or 6 columns: -- #READINFO (AllReads): total number of reads mapped -- #READINFO (UsableReads): total number of de-duplicated mapped reads -- #READINFO (GenomicReads): total number of uniquely mapped genomic reads -- #READINFO (RepFamilyReads): total number of repetitive family mapped reads -- total_or_element: this is either TOTAL or ELEMENT - - Total: family (you can create this by summing all of the elements within this family) - - total - - family name - - number of reads - - reads per million - - Element: element - - element - - family name (eg. family1). Can have multiple pipe-delimited families - - number of reads - - reads per million - - family1||transcript1|transcript2|transcript3 (transcript 1/2/3 are all members of family1) - - gene name of the transcripts -.reparsed file: tabbed file with similar structure to .parsed, but without the (Allreads) total. -.nopipes.tsv: tabbed file containing only non-ambiguously mapped families with the following columns: -- element: the repeat element -- IP_read_num: number of reads mapping to the element in IP -- IP_clip_rpr: number of mapped reads / number of total usable mapped reads -- Input_read_num: number of reads mapping to the element in Input. +1 pseudocount if no input reads were mapped. -- Input_clip_rpr: number of mapped input reads / number of total usable mapped input reads -- Fold_enrichment: IP_clip_rpr / Input_clip_rpr -- Information_content: IP_clip_rpr * log2(IP_clip_rpr / Input_clip_rpr) -.withpipes.tsv: tabbed file containing ambiguously and non-ambiguously mapped families - -# Notes: -- Elements (second column) that contain pipes ```|``` are multifamily mapped -which means that the read is ambiguously mapped to each element. These can generally -be excluded from downstream analysis. -- You can sort by information content and log2 fold enrichment. -Information content is calculated as: log2(clip_rpr/input_rpr). -- High fold changes and information content (no guidelines yet to this cutoff) -typically indicate elements enriched in the experiment. - -# References: -- Van Nostrand, E.L., Pratt, G.A., Yee, B.A. et al. Principles of RNA processing from analysis of enhanced CLIP maps for 150 RNA binding proteins. Genome Biol 21, 90 (2020). https://doi.org/10.1186/s13059-020-01982-9 +### Workflow validation run +```bash +HOME=$(pwd) ./.conda-env/bin/snakemake -j1 -p -F all +``` + +## Regenerating expected artifacts +If you intentionally update core Python logic, regenerate mini expected artifacts and commit them: +```bash +python3 scripts/generate_python_expected_artifacts.py +``` +This refreshes: +- `tests/fixtures/mini/expected/split.expected.manifest.tsv` +- `tests/fixtures/mini/expected/merged.expected.parsed` diff --git a/Snakefile b/Snakefile new file mode 100644 index 0000000..fc73a14 --- /dev/null +++ b/Snakefile @@ -0,0 +1,9 @@ +configfile: 'config/config.yaml' + +from workflow.config_adapter import merge_cwl_job_yaml_into_config + + +if config.get('cwl_input_yaml'): + config = merge_cwl_job_yaml_into_config(config, config['cwl_input_yaml']) + +include: 'workflow/rules/se_foundation.smk' diff --git a/bin/perl/._reparse_samfile_updatedchrM_fixmultenstsort_SE.pl b/bin/perl/._reparse_samfile_updatedchrM_fixmultenstsort_SE.pl deleted file mode 100644 index 236d0f9..0000000 Binary files a/bin/perl/._reparse_samfile_updatedchrM_fixmultenstsort_SE.pl and /dev/null differ diff --git a/bin/perl/RepElement_pipeline_1dataset.pl b/bin/perl/RepElement_pipeline_1dataset.pl deleted file mode 100755 index 633d91a..0000000 --- a/bin/perl/RepElement_pipeline_1dataset.pl +++ /dev/null @@ -1,173 +0,0 @@ -use warnings; -use strict; - -my $species = "hg38"; -## note - this pipeline right now is specific for hg38 - -my $keep_temp_files_flag = 0; - -## -unless ($ARGV[4]) { - print STDERR "usage: perl RepElement_pipeline_1dataset.pl working_directory fastq1_R1:fastq1_R2::fastq2_R1:fastq2_R2::(etc) genomic_mapping_bam output_prefix se_or_pe_flag\n"; - print STDERR "e.g.\n"; - print STDERR "PE example: perl RepElement_pipeline_1dataset.pl /home/elvannostrand/scratch/20201211_PEtest_204/ /projects/ps-yeolab4/software/eclip/0.7.0/tests/other_examples/204_RBFOX2_GRCh38/results/204_RBFOX2.rep1_clip.A01.r1.fqTrTr.sorted.fq.gz:/projects/ps-yeolab4/software/eclip/0.7.0/tests/other_examples/204_RBFOX2_GRCh38/results/204_RBFOX2.rep1_clip.A01.r2.fqTrTr.sorted.fq.gz::/projects/ps-yeolab4/software/eclip/0.7.0/tests/other_examples/204_RBFOX2_GRCh38/results/204_RBFOX2.rep1_clip.B06.r1.fqTrTr.sorted.fq.gz:/projects/ps-yeolab4/software/eclip/0.7.0/tests/other_examples/204_RBFOX2_GRCh38/results/204_RBFOX2.rep1_clip.B06.r2.fqTrTr.sorted.fq.gz /projects/ps-yeolab4/software/eclip/0.7.0/tests/other_examples/204_RBFOX2_GRCh38/results/204_RBFOX2.rep1_clip.A01.r1.fq.genome-mappedSo.bam::/projects/ps-yeolab4/software/eclip/0.7.0/tests/other_examples/204_RBFOX2_GRCh38/results/204_RBFOX2.rep1_clip.B06.r1.fq.genome-mappedSo.bam 204_01 PE\n"; - print STDERR "SE example: perl RepElement_pipeline_1dataset.pl /home/elvannostrand/scratch/20201211_SEtest/ /projects/ps-yeolab4/software/eclip/0.7.0/tests/other_examples/INV_B_singleNode/results/INV_B.IP.umi.r1.fqTrTr.sorted.fq.gz /projects/ps-yeolab4/software/eclip/0.7.0/tests/other_examples/INV_B_singleNode/results/INV_B.IP.umi.r1.fq.genome-mappedSoSo.bam RBFOX_seCLIP SE\n"; - exit; -} - -my $working_directory = $ARGV[0]; -my $fastqs_input = $ARGV[1]; -my $standard_analysis_bams_input = $ARGV[2]; -my $dataset_label = $ARGV[3]; -my $se_or_pe_flag = $ARGV[4]; -my @fastq_demux_pairs = split(/\:\:/,$fastqs_input); -my @standard_analysis_bams_split = split(/\:\:/,$standard_analysis_bams_input); - -my $rmdup_output_merged_sam = $working_directory.$dataset_label.".combined_w_uniquemap.rmDup.sam"; -my $merged_rmdup = $working_directory.$dataset_label."_merged.rmDup.sam"; -my $parsed_out = $rmdup_output_merged_sam.".parsed"; - -if (-e $parsed_out) { - print STDERR "skipping $dataset_label , $parsed_out already exists\n"; - next; -} - -if (-d $working_directory) { -} else { - system("mkdir $working_directory"); -} -$working_directory =~ s/\/$//; -$working_directory .= "/"; - -my $sh_out = "RM_".$dataset_label.".sh"; -open(SH,">$sh_out"); -print SH "\#\!\/bin\/sh\n"; -print SH "#PBS -N ".$sh_out."\n"; -print SH "#PBS -o ".$sh_out.".out\n"; -print SH "#PBS -e ".$sh_out.".err\n"; -print SH "#PBS -V\n"; -print SH "#PBS -l walltime=8:00:00\n"; -print SH "#PBS -l nodes=1:ppn=4\n"; -print SH "#PBS -A yeo-group\n"; -print SH "#PBS -q condo\n"; -# print SH "#PBS -q home-yeo\n"; -print SH "#PBS -M elvannostrand\@ucsd.edu\n"; -print SH "#PBS -m a\n"; -print SH "cd /home/elvannostrand/data/clip/CLIPseq_analysis/scripts/inline_processing/hg38/\n"; - - -my @unmerged_sams; -my @unmerged_parsed; - -for my $fastq (@fastq_demux_pairs) { - my @fastqpairs = split(/\:/,$fastq); - if (scalar(@fastqpairs) > 2) { - print STDERR "error - this should be either 1 fastq (single-end eCLIP) or 2 (R1 and R2 for original paired-end eCLIP), but got ".scalar(@fastqpairs)." instead\n$fastq\n"; - exit; - } - - foreach (@fastqpairs) { - my @fastq_split = split(/\//,$_); - my $fastq_short = $fastq_split[$#fastq_split]; - my $newpath = $working_directory.$fastq_short; - - print SH "cp $_ $newpath\n"; - print SH "gunzip $newpath\n"; - - $fastq_short =~ s/\.gz$//; - $_ = $working_directory.$fastq_short; - } - my $fastq1 = $fastqpairs[0]; - my $fastq2; - $fastq2 = $fastqpairs[1] if ($se_or_pe_flag eq "PE"); - - my @fastq_fi1_split = split(/\//,$fastq1); - my $fastq_fi1_short = $fastq_fi1_split[$#fastq_fi1_split]; - my $merged_bam_file_standardanalysis = shift(@standard_analysis_bams_split); - my @merged_bam_file_split = split(/\//,$merged_bam_file_standardanalysis); - my $merged_bam_file_standardanalysis_short = $merged_bam_file_split[$#merged_bam_file_split]; - - print SH "perl split_bam_to_subfiles_SEorPE.pl $merged_bam_file_standardanalysis $working_directory $se_or_pe_flag\n"; - -# my $bowtie_db = "/home/elvannostrand/data/clip/CLIPseq_analysis/RNA_type_analysis/MASTER_filelist.wrepbaseandtRNA.fa.fixed.fa.UpdatedSimpleRepeat"; - #hg19 - my $bowtie_db = "/home/elvannostrand/data/clip/CLIPseq_analysis/RNA_type_analysis/MASTER_filelist.wrepbaseandtRNA.fa.fixed.fa.UpdatedSimpleRepeat.20190424"; - if ($species eq "hg38") { - $bowtie_db = "/home/elvannostrand/data/clip/CLIPseq_analysis/RNA_type_analysis/hg38/MASTER_FILELIST.20201203.wrepbaseandtRNA.fa.fixed.fa.UpdatedSimpleRepeat"; - } -# my $bowtie_db = "/home/elvannostrand/data/clip/CLIPseq_analysis/RNA_type_analysis/MASTER_filelist.wrepbaseandtRNA.fa.fixed"; - my @bowtie_db_split = split(/\//,$bowtie_db); - my $bowtie_db_short = $bowtie_db_split[$#bowtie_db_split]; - - my $sam_output = $working_directory.$fastq_fi1_short.".mapped_vs_".$bowtie_db_short.".sam"; - my $sam_output_short = $fastq_fi1_short.".mapped_vs_".$bowtie_db_short.".sam"; - - if ($se_or_pe_flag eq "PE") { - print SH "perl parse_bowtie2_output_realtime_includemultifamily_PE.pl $fastq1 $fastq2 $bowtie_db $working_directory $sam_output $species\n"; - } else { - print SH "perl parse_bowtie2_output_realtime_includemultifamily_SE.pl $fastq1 $bowtie_db $working_directory $sam_output $species\n"; - } - print SH "rm $fastq1\n" unless ($keep_temp_files_flag==1); - if ($se_or_pe_flag eq "PE") { - print SH "rm $fastq2\n" unless ($keep_temp_files_flag==1); - } - - # now split sam file into files separated by first 2 bases of barcode - print SH "perl split_bam_to_subfiles_SEorPE.pl $sam_output $working_directory $se_or_pe_flag\n"; - - print SH "rm $sam_output\n" unless ($keep_temp_files_flag==1); - - my $predup_output = $working_directory.$sam_output_short.".combined_w_uniquemap.prermDup.sam"; - - my $rmdup_output = $working_directory.$sam_output_short.".combined_w_uniquemap.rmDup.sam"; - print SH "rm $predup_output\n"; - print SH "rm $rmdup_output\n"; -# print SH "rm $rmdup_output\n" if (-e $rmdup_output); - - for my $b1 ("A","C","G","T","N") { - for my $b2 ("A","C","G","T","N") { - my $repmapping_fi = $working_directory.$sam_output_short.".".$b1.$b2.".tmp"; - my $genomemapping_fi = $working_directory.$merged_bam_file_standardanalysis_short.".".$b1.$b2.".tmp"; - - print SH "perl duplicate_removal_inline_paired.count_region_other_reads_masksnRNAs_andreparse_SEandPE_20201210_simple.pl $repmapping_fi $genomemapping_fi $se_or_pe_flag\n"; - my $tmp_output_fi = $working_directory.$sam_output_short.".".$b1.$b2.".tmp".".combined_w_uniquemap.rmDup.sam"; - my $tmp_parsed_fi = $working_directory.$sam_output_short.".".$b1.$b2.".tmp".".combined_w_uniquemap.rmDup.sam.parsed_v2.20201210.txt"; - my $tmp_predup_fi = $working_directory.$sam_output_short.".".$b1.$b2.".tmp".".combined_w_uniquemap.prermDup.sam"; - - print SH "cat $tmp_output_fi >> $rmdup_output\n"; - print SH "cat $tmp_predup_fi >> $predup_output\n"; - - print SH "rm $repmapping_fi\n" unless ($keep_temp_files_flag==1); - print SH "rm $genomemapping_fi\n" unless ($keep_temp_files_flag==1); - print SH "rm $tmp_output_fi\n" unless ($keep_temp_files_flag==1); - print SH "rm $tmp_predup_fi\n" unless ($keep_temp_files_flag==1); - push @unmerged_parsed,$tmp_parsed_fi; - } - } - - print SH "gzip $predup_output\n"; - push @unmerged_sams,$rmdup_output; -} - -print SH "perl merge_multiple_parsed_files.simplified_20191022.pl $parsed_out ".join(" ",@unmerged_parsed)."\n"; -for my $unmerged_parsed_fi (@unmerged_parsed) { - print SH "rm $unmerged_parsed_fi\n" unless ($keep_temp_files_flag==1); -} - -print SH "rm $rmdup_output_merged_sam\n"; -my $merge_command = "cp $unmerged_sams[0] $rmdup_output_merged_sam\n"; -for (my $i=1;$i families) - - -my %revstrand; -$revstrand{"+"} = "-"; -$revstrand{"-"} = "+"; - -my $genome_hashing_value = 1000; -my %enst2type; -my %enst2ensg; -my %ensg2name; -my %ensg2type; - - -my $gencode_gtf_file = $ARGV[3]; -my $gencode_tablebrowser_file = $ARGV[4]; -my $repmask_bed_fi = $ARGV[5]; -my $filelist_file = $ARGV[6]; - -my %gencode_features; -&read_gencode_gtf($gencode_gtf_file); -&read_gencode($gencode_tablebrowser_file); -# this flag refers to how to count uniquely genomic mapped reads - 5' end only? -my $region_5primeend_only_flag = 1; - - - - -my %enst2gene; -my %convert_enst2type; -my %peaks; -&read_peakfi($repmask_bed_fi); - -my %convert_enst2priorityN; -my %convert_enst2priority; - -&read_in_filelists($filelist_file); - -#&read_in_filelists($filelist_file2); - -#my %mirbase_features; - -#&read_mirbase($mirbase_fi); -my %convert_strand = ("+" => "-", "-" => "+"); - - -my $unique_count = 0; -my ($all_count,$duplicate_count,$unique_count_old,$unique_genomic_count,$unique_repfamily_count) = (0,0,0,0,0); - -# takes in sam file, should be blah.bc.tmp file -my $repfamily_sam = $ARGV[0]; -my @repfamily_sam_split = split(/\//,$repfamily_sam); -my $repfamily_sam_short = $repfamily_sam_split[$#repfamily_sam_split]; - -my $gabe_rmRep_sam = $ARGV[1]; -my $eCLIP_read_type_flag = $ARGV[2]; - - -my %read_hash; -my $output_fi = $repfamily_sam_short.".combined_w_uniquemap.rmDup.sam"; -open(OUT,">$output_fi"); - -my ($total_unique_mapped_read_num,$rep_family_reads,$unique_genomic) = (0,0,0); - -my %count; -my %count_enst; -my $pre_rmdup_fi = $repfamily_sam_short.".combined_w_uniquemap.prermDup.sam"; -open(PREDUP,">$pre_rmdup_fi"); - - -if ($eCLIP_read_type_flag eq "SE") { - print STDERR "running in single-end mode\n"; - &read_rep_family_se($repfamily_sam); - &read_unique_mapped_se($gabe_rmRep_sam); - &run_pcr_duplicate_removal_se(); -} elsif ($eCLIP_read_type_flag eq "PE") { - print STDERR "running in paired-end mode\n"; - &read_rep_family_pe($repfamily_sam); - &read_unique_mapped_pe($gabe_rmRep_sam); - &run_pcr_duplicate_removal_pe(); -} - - - -close(PREDUP); -close(OUT); -my $count_out = $output_fi.".parsed_v2.20201210.txt"; -open(COUNT,">$count_out"); -print COUNT "#READINFO\tAll reads:\t$all_count\tPCR duplicates removed:\t$duplicate_count\tUsable Remaining:\t$unique_count_old\tUsable from genomic mapping:\t$unique_genomic_count\tUsable from family mapping:\t$unique_repfamily_count\n"; -#print COUNT "#READINFO\tAll reads:\t$all_count\nPCR duplicates removed:\t$duplicate_count\nUsable Remainig:\t$unique_count\nUsable from genomic mapping:\t$unique_genomic_count\nUsable from family mapping:\t$unique_repfamily_count\n"; - -#&count_output(); -print COUNT "#READINFO\tUsableReads\t".$total_unique_mapped_read_num."\n"; -if ($total_unique_mapped_read_num > 0) { - print COUNT "#READINFO\tGenomicReads\t".$unique_genomic."\t".sprintf("%.5f",$unique_genomic/$total_unique_mapped_read_num)."\n"; - print COUNT "#READINFO\tRepFamilyReads\t".$rep_family_reads."\t".sprintf("%.5f",$rep_family_reads/$total_unique_mapped_read_num)."\n"; -} else { - print COUNT "#READINFO\tGenomicReads\t".$unique_genomic."\t0\n"; - print COUNT "#READINFO\tRepFamilyReads\t".$rep_family_reads."\t0\n"; -} - -my @sorted_total = sort {$count{$b} <=> $count{$a}} keys %count; -for my $k (@sorted_total) { - print COUNT "TOTAL\t$k\t$count{$k}\t".sprintf("%.5f",$count{$k} * 1000000 / $total_unique_mapped_read_num)."\n"; -} - -my @sorted = sort {$count_enst{$b} <=> $count_enst{$a}} keys %count_enst; -for my $s (@sorted) { - my ($ensttype,$multensts) = split(/\|\|/,$s); - my @multensts_split = split(/\|/,$multensts); - my @genes_final; - my $type = $ensttype; - for my $multensts_group (@multensts_split) { - my @gids = split(/\;\;/,$multensts_group); - - my @genes_short; - for my $gid (@gids) { - if (exists $enst2gene{$gid}) { - push @genes_short,$enst2gene{$gid}; - } else { - push @genes_short,$gid; - } - } - push @genes_final,join(";;",@genes_short); - } - - my ($ensg_primary,$readnum,$rpm,$enst_all,$ensg_all) = ($type,$count_enst{$s},sprintf("%.5f",$count_enst{$s} * 1000000 / $total_unique_mapped_read_num),$s,join("|",@genes_final)); - - - print COUNT "ELEMENT\t".$ensg_primary."\t".$readnum."\t".$rpm."\t".$enst_all."\t".$ensg_all."\n"; -} - - - -close(COUNT); - - - - -my $out_done = $count_out.".done"; -open(DONE,">$out_done"); -print DONE "jobs done\n"; -close(DONE); - - -# now do PCR duplicate removal - -sub run_pcr_duplicate_removal_se { - my %fragment_hash; - -#changed 20171019 to sort reads to be stable across perl versions - my @sorted_reads = sort {$a cmp $b} keys %read_hash; - for my $r1name (@sorted_reads) { -# for my $r1name (keys %read_hash) { - my $r1 = $read_hash{$r1name}{R1}; - - my @tmp_r1 = split(/\t/,$r1); - my ($r1name,$r1bc) = split(/\s+/,$tmp_r1[0]); - - my $r1sam_flag = $tmp_r1[1]; - next if ($r1sam_flag == 4); -# next if ($r1sam_flag == 77 || $r1sam_flag == 141); - - my $frag_strand; -### This section is for only properly paired reads -# if ($r1sam_flag == 99 || $r1sam_flag == 355) { -# $frag_strand = "-"; -# } elsif ($r1sam_flag == 83 || $r1sam_flag == 339) { -# $frag_strand = "+"; -# } elsif ($r1sam_flag == 147 || $r1sam_flag == 403) { -# $frag_strand = "-"; -# } elsif ($r1sam_flag == 163 || $r1sam_flag == 419) { -# $frag_strand = "+"; -# } else { -# next; -# print STDERR "R1 strand error $r1sam_flag\n"; -# } - - if ($r1sam_flag == 16 || $r1sam_flag == 272) { - $frag_strand = "-"; - } elsif ($r1sam_flag eq "0" || $r1sam_flag == 256) { - $frag_strand = "+"; - } else { - next; - print STDERR "R1 strand error $r1sam_flag\n"; - } - - - my @read_name = split(/\_/,$tmp_r1[0]); - my $randommer = pop(@read_name); - - my $r1_cigar = $tmp_r1[5]; - - # 165 = R2 unmapped, R1 rev strand -- frag on fwd strand - - - my $r1_chr = $tmp_r1[2]; - my $r1_start = $tmp_r1[3]; - - my $mismatch_flags = $tmp_r1[14]; - - my %read_regions; - @{$read_regions{"R1"}} = &parse_cigar_string($r1_start,$r1_cigar,$r1_chr,$frag_strand); - - my $hashing_value; - if ($frag_strand eq "+") { - my $r1_firstregion = $read_regions{"R1"}[scalar(@{$read_regions{"R1"}})-1]; - my ($r1_firstchr,$r1_firststr,$r1_firstpos) = split(/\:/,$r1_firstregion); - my ($r1_firststart,$r1_firststop) = split(/\-/,$r1_firstpos); - $hashing_value = $r1_firstchr.":".$r1_firststr.":".$r1_firststart."\t".$r1_firstchr.":".$r1_firststr.":".$r1_firststop; - } elsif ($frag_strand eq "-") { - my $r1_firstregion = $read_regions{"R1"}[0]; - my ($r1_firstchr,$r1_firststr,$r1_firstpos) = split(/\:/,$r1_firstregion); - my ($r1_firststart,$r1_firststop) = split(/\-/,$r1_firstpos); - $hashing_value = $r1_firstchr.":".$r1_firststr.":".$r1_firststart."\t".$r1_firstchr.":".$r1_firststr.":".$r1_firststop; - } else { - print STDERR "strand error $frag_strand $r1\n"; - } - - my $full_frag_position = join("|",@{$read_regions{"R1"}}); - - my $debug_flag = 0; - - my $hashing_value_toprint = $r1_chr."|".$frag_strand."::".$hashing_value.":".$randommer; - print PREDUP "$r1\t".$hashing_value_toprint."\n"; - - my %full_fragment; - $all_count++; - if (exists $fragment_hash{$r1_chr."|".$frag_strand}{$hashing_value.":".$randommer}) { - $duplicate_count++; - delete($read_hash{$r1name}); - next; - } else { - $fragment_hash{$r1_chr."|".$frag_strand}{$hashing_value.":".$randommer} = 1; - - if (exists $read_hash{$r1name}{file1flag} && $read_hash{$r1name}{file1flag} == 1) { - # read is coming from familiy mapping - print OUT "".$r1."\tRepFamily\t$r1_chr\n"; - push @tmp_r1,"RepFamily"; -# push @tmp_r1,$r1_chr; - push @tmp_r1,$read_hash{$r1name}{ensttype}."||".$read_hash{$r1name}{mult_ensts}; - - $unique_repfamily_count++; - } elsif (exists $read_hash{$r1name}{file2flag} && $read_hash{$r1name}{file2flag} == 1) { - # read is coming from genomic mapping - print OUT "".$r1."\tUniqueGenomic\t".$read_hash{$r1name}{ensttype}."||".$read_hash{$r1name}{mult_ensts}."\n"; - push @tmp_r1,"UniqueGenomic"; - push @tmp_r1,$read_hash{$r1name}{ensttype}."||".$read_hash{$r1name}{mult_ensts}; - - $unique_genomic_count++; - } else { - print STDERR "this shouldn't be hit - fi1flag $read_hash{$r1name}{file1flag} fi2flag $read_hash{$r1name}{file2flag} $r1name\n"; - } - } - - $unique_count_old++; - - - my $repmap_info = pop(@tmp_r1); - my $type = pop(@tmp_r1); - my $mult_transcripts = pop(@tmp_r1); - - my ($ensttype,$mult_ensts) = split(/\|\|/,$repmap_info); - - if ($type eq "RepFamily") { - $rep_family_reads++; - } elsif ($type eq "UniqueGenomic") { - if ($r1_chr eq "chrM") { - $ensttype = "chrM_unique_".$frag_strand."strand"; - } - $unique_genomic++; - } else { - print STDERR "error - this shouldn't happen $r1\n"; - } - - if ($ensttype =~ /Simple\_repeat/) { - $ensttype = "Simple_repeat"; - } - - $count{$ensttype}++; - $count_enst{$ensttype."||".$mult_ensts}++; - $total_unique_mapped_read_num++ - - } -} - -sub run_pcr_duplicate_removal_pe { - my %fragment_hash; - -#changed 20171019 to sort reads to be stable across perl versions - my @sorted_reads = sort {$a cmp $b} keys %read_hash; - for my $r1name (@sorted_reads) { -# for my $r1name (keys %read_hash) { - my $r1 = $read_hash{$r1name}{R1}; - my $r2 = $read_hash{$r1name}{R2}; - - my @tmp_r1 = split(/\t/,$r1); - my @tmp_r2 = split(/\t/,$r2); - my ($r1name,$r1bc) = split(/\s+/,$tmp_r1[0]); - my ($r2name,$r2bc) = split(/\s+/,$tmp_r2[0]); - - my $r1sam_flag = $tmp_r1[1]; - my $r2sam_flag = $tmp_r2[1]; - unless ($r1sam_flag) { - print STDERR "error $r1 $r2\n"; - } - next if ($r1sam_flag == 77 || $r1sam_flag == 141); - - my $frag_strand; -### This section is for only properly paired reads - if ($r1sam_flag == 99 || $r1sam_flag == 355) { - $frag_strand = "-"; - } elsif ($r1sam_flag == 83 || $r1sam_flag == 339) { - $frag_strand = "+"; - } elsif ($r1sam_flag == 147 || $r1sam_flag == 403) { - $frag_strand = "-"; - } elsif ($r1sam_flag == 163 || $r1sam_flag == 419) { - $frag_strand = "+"; - } else { - next; - print STDERR "R1 strand error $r1sam_flag\n"; - } - - my @read_name = split(/\:/,$tmp_r1[0]); - my $randommer = $read_name[0]; - - my $r1_cigar = $tmp_r1[5]; - my $r2_cigar = $tmp_r2[5]; - - # 165 = R2 unmapped, R1 rev strand -- frag on fwd strand - - - my $r1_chr = $tmp_r1[2]; - my $r1_start = $tmp_r1[3]; - my $r2_chr = $tmp_r2[2]; - my $r2_start = $tmp_r2[3]; - - my $mismatch_flags = $tmp_r1[14]; - - my %read_regions; - @{$read_regions{"R1"}} = &parse_cigar_string($r1_start,$r1_cigar,$r1_chr,$frag_strand); - @{$read_regions{"R2"}} = &parse_cigar_string($r2_start,$r2_cigar,$r2_chr,$frag_strand); - - my $hashing_value; - if ($frag_strand eq "+") { - my $r1_firstregion = $read_regions{"R1"}[scalar(@{$read_regions{"R1"}})-1]; - my ($r1_firstchr,$r1_firststr,$r1_firstpos) = split(/\:/,$r1_firstregion); - my ($r1_firststart,$r1_firststop) = split(/\-/,$r1_firstpos); - my $r2_firstregion = $read_regions{"R2"}[0]; - my ($r2_firstchr,$r2_firststr,$r2_firstpos) = split(/\:/,$r2_firstregion); - my ($r2_firststart,$r2_firststop) = split(/\-/,$r2_firstpos); - $hashing_value = $r2_firstchr.":".$r2_firststr.":".$r2_firststart."\t".$r1_firstchr.":".$r1_firststr.":".$r1_firststop; - } elsif ($frag_strand eq "-") { - my $r1_firstregion = $read_regions{"R1"}[0]; - my ($r1_firstchr,$r1_firststr,$r1_firstpos) = split(/\:/,$r1_firstregion); - my ($r1_firststart,$r1_firststop) = split(/\-/,$r1_firstpos); - my $r2_firstregion = $read_regions{"R2"}[scalar(@{$read_regions{"R2"}})-1]; - my ($r2_firstchr,$r2_firststr,$r2_firstpos) = split(/\:/,$r2_firstregion); - my ($r2_firststart,$r2_firststop) = split(/\-/,$r2_firstpos); - $hashing_value = $r1_firstchr.":".$r1_firststr.":".$r1_firststart."\t".$r2_firstchr.":".$r2_firststr.":".$r2_firststop; - - } else { - print STDERR "strand error $frag_strand $r1 $r2\n"; - } - my $full_frag_position = join("|",@{$read_regions{"R1"}})."\t".join("|",@{$read_regions{"R2"}}); - my $debug_flag = 0; - - my $hashing_value_toprint = $r1_chr."|".$frag_strand."::".$hashing_value.":".$randommer; - print PREDUP "$r1\t".$hashing_value_toprint."\n$r2\t".$hashing_value_toprint."\n"; - - my %full_fragment; - $all_count++; - if (exists $fragment_hash{$r1_chr."|".$frag_strand}{$hashing_value.":".$randommer}) { - $duplicate_count++; - delete($read_hash{$r1name}); - next; - } else { - $fragment_hash{$r1_chr."|".$frag_strand}{$hashing_value.":".$randommer} = 1; - - if (exists $read_hash{$r1name}{file1flag} && $read_hash{$r1name}{file1flag} == 1) { - # read is coming from familiy mapping - print OUT "".$r1."\tRepFamily\t$r1_chr\n".$r2."\tRepFamily\t$r1_chr\n"; - push @tmp_r1,"RepFamily"; -# push @tmp_r1,$r1_chr; - push @tmp_r1,$read_hash{$r1name}{ensttype}."||".$read_hash{$r1name}{mult_ensts}; - - $unique_repfamily_count++; - } elsif (exists $read_hash{$r1name}{file2flag} && $read_hash{$r1name}{file2flag} == 1) { - # read is coming from familiy mapping - print OUT "".$r1."\tUniqueGenomic\t".$read_hash{$r1name}{ensttype}."||".$read_hash{$r1name}{mult_ensts}."\n".$r2."\tUniqueGenomic\t".$read_hash{$r1name}{ensttype}."||".$read_hash{$r1name}{mult_ensts}."\n"; - push @tmp_r1,"UniqueGenomic"; - push @tmp_r1,$read_hash{$r1name}{ensttype}."||".$read_hash{$r1name}{mult_ensts}; - - $unique_genomic_count++; - } else { - print STDERR "this shouldn't be hit - fi1flag $read_hash{$r1name}{file1flag} fi2flag $read_hash{$r1name}{file2flag} $r1name\n"; - } - } - - $unique_count_old++; - my $repmap_info = pop(@tmp_r1); - my $type = pop(@tmp_r1); - my $mult_transcripts = pop(@tmp_r1); - - my ($ensttype,$mult_ensts) = split(/\|\|/,$repmap_info); - - if ($type eq "RepFamily") { - $rep_family_reads++; - } elsif ($type eq "UniqueGenomic") { - if ($r1_chr eq "chrM") { - $ensttype = "chrM_unique_".$frag_strand."strand"; - } - $unique_genomic++; - } else { - print STDERR "error - this shouldn't happen $r1\n"; - } - - if ($ensttype =~ /Simple\_repeat/) { - $ensttype = "Simple_repeat"; - } - - $count{$ensttype}++; - $count_enst{$ensttype."||".$mult_ensts}++; - $total_unique_mapped_read_num++ - - } -} - - - - -sub read_unique_mapped_pe { - my %deleted; - my $fi2_count=0; - my %read_hash_bam2; - my $gabe_rmRep_sam_file = shift; -# print STDERR "opening $gabe_rmRep_sam_file\n"; - if ($gabe_rmRep_sam_file =~ /\.sam$/ || $gabe_rmRep_sam_file =~ /\.tmp$/) { - open(B,$gabe_rmRep_sam_file) || die "no $gabe_rmRep_sam_file\n"; - } elsif ($gabe_rmRep_sam_file =~/\.bam$/) { - open(B,"samtools view -h $gabe_rmRep_sam_file |") || die "no $gabe_rmRep_sam_file\n"; - } else { - die "couldn't figure out format of $gabe_rmRep_sam_file\n"; - } - while () { - my $r1 = $_; - next if ($r1 =~ /^\@/); - $fi2_count++; - print STDERR "read $fi2_count\n" if ($fi2_count % 100000 == 0); - my $r2 = ; - chomp($r1); - chomp($r2); - - my @tmp_r1 = split(/\t/,$r1); - my @tmp_r2 = split(/\t/,$r2); - my ($r1name,$r1bc) = split(/\s+/,$tmp_r1[0]); - my ($r2name,$r2bc) = split(/\s+/,$tmp_r2[0]); - - unless ($tmp_r1[0] eq $tmp_r2[0]) { - print STDERR "paired end mismatch error: $tmp_r1[0] $tmp_r2[0]\n"; - } - - my $r2sam_flag = $tmp_r2[1]; - my $r1sam_flag = $tmp_r1[1]; - next if ($r1sam_flag == 77 || $r1sam_flag == 141); - - - # 77 = R1, unmapped - # 141 = R2, unmapped - - # 99 = R1, mapped, fwd strand --- frag on rev strand - # 101 = R1 unmapped, R2 mapped rev strand -- frag on rev strand - # 73 = R1, mapped, fwd strand --- frag on rev strand - # 147 = R2, mapped, rev strand -- frag on rev strand - # 153 = R2 mapped (R1 unmapped), rev strand -- frag on rev strand - # 133 = R2 unmapped, R1 mapped fwd strand -- frag on rev strand - - # 83 = R1, mapped, rev strand --- frag on fwd strand - # 69 = R1 unmapped, R2 mapped fwd strand -- frag on fwd strand - # 89 = R1 mapped rev strand, R2 unmapped -- frag on fwd strand - # 163 = R2, mapped, fwd strand -- frag on fwd strand - # 137 = R2 mapped (R1 unmapped), fwd strand -- frag on fwd strand - # 165 = R2 unmapped, R1 rev strand -- frag on fwd strand - - - my $frag_strand; -### This section is for only properly paired reads - if ($r1sam_flag == 99) { - $frag_strand = "-"; - } elsif ($r1sam_flag == 83) { - $frag_strand = "+"; - } elsif ($r1sam_flag == 147) { - $frag_strand = "-"; - ($r1,$r2) = ($r2,$r1); - @tmp_r1 = split(/\t/,$r1); - @tmp_r2 = split(/\t/,$r2); - - } elsif ($r1sam_flag == 163) { - ($r1,$r2) =($r2,$r1); - @tmp_r1 = split(/\t/,$r1); - @tmp_r2 = split(/\t/,$r2); - $frag_strand = "+"; - } else { -# print STDERR "R1 strand error $r1sam_flag\n"; - next; - } - my @read_name = split(/\:/,$tmp_r1[0]); - my $randommer = $read_name[0]; - - my $r1_cigar = $tmp_r1[5]; - my $r2_cigar = $tmp_r2[5]; - - my $r1_chr = $tmp_r1[2]; - my $r1_start = $tmp_r1[3]; - - my $r2_chr = $tmp_r2[2]; - my $r2_start = $tmp_r2[3]; - - - my $flags_r1 = join("\t",@tmp_r1[11..$#tmp_r1]); - my $flags_r2 = join("\t",@tmp_r2[11..$#tmp_r2]); - - my $r1_mismatch; - my $r2_mismatch; - if ($flags_r1 =~ /MD\:Z\:(\S+)\s/ || $flags_r1 =~ /MD\:Z\:(\S+?)$/) { - $r1_mismatch = $1; - } - if ($flags_r2 =~ /MD\:Z\:(\S+)\s/ || $flags_r2 =~ /MD\:Z\:(\S+?)$/) { - $r2_mismatch = $1; - } - - my $r1_phred = $tmp_r1[10]; - my $r2_phred = $tmp_r2[10]; - my $r1_seq = $tmp_r1[9]; - my $r2_seq = $tmp_r2[9]; - my $r1_mmscore = &get_alignment_score($r1_mismatch,$r1_cigar,$r1_phred,$r1_seq); - my $r2_mmscore = &get_alignment_score($r2_mismatch,$r2_cigar,$r2_phred,$r2_seq); - my $total_mmscore = $r1_mmscore+$r2_mmscore; - - if (exists $read_hash{$r1name}{file1flag}) { - unless (exists $read_hash{$r1name}{rep_score}) { - print STDERR "err $r1name\n"; - } - if ($total_mmscore > $read_hash{$r1name}{rep_score} + 24) { - - # if unique mapping to genome is more than 2 mismatches per read = 2 * 2 * 6 alignment score better than to repeat element, throw out repeat element and use genome mapping -# print "deleting $r1name from rep element mapping - $total_mmscore $read_hash{$r1name}{rep_score}\n$r1\n$read_hash{$r1name}{R1}\n----\n"; - $deleted{$read_hash{$r1name}{ensttype}}++; - delete($read_hash{$r1name}) if (exists $read_hash{$r1name}); - - } else { - # repeat mapping exists and is better than unique genome mapping - throw out unique genome mapping - next; - } - } - - my @read_regions = &parse_cigar_string($r2_start,$r2_cigar,$r2_chr,$frag_strand); - my %tmp_hash; - for my $region (@read_regions) { - my ($rchr,$rstr,$rpos) = split(/\:/,$region); - my ($rstart,$rstop) = split(/\-/,$rpos); - - my $verbose_flag = 0; - - my $rx = int($rstart / $genome_hashing_value); - my $ry = int($rstop / $genome_hashing_value); - for my $ri ($rx..$ry) { - - for my $peak (@{$peaks{$rchr}{$rstr}{$ri}}) { - my ($pchr,$ppos,$pstr,$ptype) = split(/\:/,$peak); - my ($pstart,$pstop) = split(/\-/,$ppos); - -# next if ($pstart > $rstop || $pstop < $rstart); -# if ($pstop == $rstart) { -# if ($pstart == $rstop) { -# if ($pstart == $rstart) { -# print "region $region $r1 peak $peak\n"; -# } - next if ($pstart >= $rstop || $pstop <= $rstart); - $tmp_hash{$peak} = 1; - - } - } - } - - my %temp_peak_read_counts; - my %converted_types; - for my $peak (keys %tmp_hash) { - my ($pchr,$ppos,$pstr,$ptype) = split(/\:/,$peak); - $temp_peak_read_counts{$ptype."_uniquegenomic"}++; - $converted_types{$convert_enst2type{$ptype}} = 1; - print STDERR "peak $peak $ptype\n" unless (exists $convert_enst2type{$ptype}); - } - - - my @sorted_types = sort {$a cmp $b} keys %temp_peak_read_counts; - my $all_mapped_ensts = join("|",@sorted_types); - my @sorted_convertedtype = sort {$a cmp $b} keys %converted_types; - my $all_ensttypes = join("|",@sorted_convertedtype); - - if (exists $converted_types{"miRNA"} || exists $converted_types{"miRNA-proximal"}) { - my $final_feature_type; - if (exists $converted_types{"miRNA"}) { - $final_feature_type = "miRNA"; - } elsif (exists $converted_types{"miRNA-proximal"}) { - $final_feature_type = "miRNA-proximal"; - } - - $all_ensttypes = "unique_".$final_feature_type; -#20171019 - changed this to sort to be consistent across perl versions -# $all_mapped_ensts = "unique_".join("|",@sorted_keys); - } else { - - my $read2_start_position; - if ($frag_strand eq "+") { - $read2_start_position = $r2_start; - } elsif ($frag_strand eq "-") { - my $last_region = $read_regions[$#read_regions]; - my ($rchr,$rstr,$rpos) = split(/\:/,$last_region); - my ($rstart,$rstop) = split(/\-/,$rpos); - $read2_start_position = $rstop - 1; - } else { - print STDERR "error $frag_strand\n"; - } - - - unless (scalar(keys %temp_peak_read_counts) >= 1) { - my $feature_flag = 0; - my %tmp_gencode_hash; - - my $rx = int($read2_start_position / $genome_hashing_value); - for my $gencode (@{$gencode_features{$r1_chr}{$frag_strand}{$rx}}) { - my ($gencode_enst,$gencode_type,$gencode_region) = split(/\|/,$gencode); - my ($gencode_start,$gencode_stop) = split(/\-/,$gencode_region); - - next if ($read2_start_position < $gencode_start); - next if ($read2_start_position >= $gencode_stop); - my $gencode_ensg = $enst2ensg{$gencode_enst}; - $tmp_gencode_hash{$gencode_type}{$gencode_ensg}="contained"; - $feature_flag = 1; - - } - my $final_feature_type = "intergenic"; - - if ($feature_flag == 1) { - if (exists $tmp_gencode_hash{"CDS"}) { - my $feature_type_flag = &get_type_flag(\%tmp_gencode_hash,"CDS"); - $final_feature_type = "CDS"; - - } elsif (exists $tmp_gencode_hash{"3utr"} || exists $tmp_gencode_hash{"5utr"}) { - if (exists $tmp_gencode_hash{"3utr"} && exists $tmp_gencode_hash{"5utr"}) { - my $feature_type_flag = &get_type_flag(\%tmp_gencode_hash,"3utr"); - my $feature_type_flag2 = &get_type_flag(\%tmp_gencode_hash,"5utr"); - - unless ($feature_type_flag eq "contained" && $feature_type_flag2 eq "contained") { - $feature_type_flag = "partial"; - } - - $final_feature_type = "5utr_and_3utr"; - } elsif (exists $tmp_gencode_hash{"3utr"}) { - my $feature_type_flag = &get_type_flag(\%tmp_gencode_hash,"3utr"); - $final_feature_type = "3utr"; - } elsif (exists $tmp_gencode_hash{"5utr"}) { - my $feature_type_flag = &get_type_flag(\%tmp_gencode_hash,"5utr"); - $final_feature_type = "5utr"; - } else { - print STDERR "weird shouldn't hit this\n"; - } - } elsif (exists $tmp_gencode_hash{"proxintron"}) { - my $feature_type_flag = &get_type_flag(\%tmp_gencode_hash,"proxintron"); - $final_feature_type = "proxintron"; - } elsif (exists $tmp_gencode_hash{"distintron"}) { - my $feature_type_flag = &get_type_flag(\%tmp_gencode_hash,"distintron"); - $final_feature_type = "distintron"; - } elsif (exists $tmp_gencode_hash{"noncoding_exon"}) { - my $feature_type_flag = &get_type_flag(\%tmp_gencode_hash,"noncoding_exon"); - $final_feature_type = "noncoding_exon"; - } elsif (exists $tmp_gencode_hash{"noncoding_proxintron"}) { - my $feature_type_flag = &get_type_flag(\%tmp_gencode_hash,"noncoding_proxintron"); - $final_feature_type = "noncoding_proxintron"; - } elsif (exists $tmp_gencode_hash{"noncoding_distintron"}) { - my $feature_type_flag = &get_type_flag(\%tmp_gencode_hash,"noncoding_distintron"); - $final_feature_type = "noncoding_distintron"; - } elsif (exists $tmp_gencode_hash{"antisense_gencode"}) { - $final_feature_type = "antisense_gencode"; - } else { - print STDERR "weird - shouldn't hit this? elements of tmp_gencode_hash are ".join("|",keys %tmp_gencode_hash)."\n"; - } - } - -# $read_counts{$final_feature_type}++; - - $all_ensttypes = "unique_".$final_feature_type; - if (exists $tmp_gencode_hash{$final_feature_type}) { - -#20171019 - changed this to sort to be consistent across perl versions - my @sorted_keys = sort {$a cmp $b} keys %{$tmp_gencode_hash{$final_feature_type}}; - $all_mapped_ensts = "unique_".join("|",@sorted_keys); -# $all_mapped_ensts = "unique_".join("|",keys %{$tmp_gencode_hash{$final_feature_type}}); - } else { - $all_mapped_ensts = "unique_".$final_feature_type; - } - } - } -# next unless (scalar(keys %temp_peak_read_counts) >= 1); - - - -# print "keeping new mapping $all_mapped_ensts \t $all_ensttypes\n$r1\n"; - - $read_hash{$r1name}{rep_flag} = 0; - $read_hash{$r1name}{rep_score} = $total_mmscore; - $read_hash{$r1name}{file2flag} = 1; - $read_hash{$r1name}{R1} = $r1; - $read_hash{$r1name}{R2} = $r2; - $read_hash{$r1name}{mult_ensts} = $all_mapped_ensts; - $read_hash{$r1name}{ensttype} = $all_ensttypes; - - } - close(B); -} - - - - - -sub read_unique_mapped_se { - my %deleted; - my $fi2_count=0; - my %read_hash_bam2; - my $gabe_rmRep_sam_file = shift; -# print STDERR "opening $gabe_rmRep_sam_file\n"; - if ($gabe_rmRep_sam_file =~ /\.sam$/ || $gabe_rmRep_sam_file =~ /\.tmp$/) { - open(B,$gabe_rmRep_sam_file) || die "no $gabe_rmRep_sam_file\n"; - } elsif ($gabe_rmRep_sam_file =~/\.bam$/) { - open(B,"samtools view -h $gabe_rmRep_sam_file |") || die "no $gabe_rmRep_sam_file\n"; - } else { - die "couldn't figure out format of $gabe_rmRep_sam_file\n"; - } - while () { - my $r1 = $_; - next if ($r1 =~ /^\@/); - $fi2_count++; - print STDERR "read $fi2_count\n" if ($fi2_count % 100000 == 0); - chomp($r1); - - my @tmp_r1 = split(/\t/,$r1); - my ($r1name,$r1bc) = split(/\s+/,$tmp_r1[0]); - - my $r1sam_flag = $tmp_r1[1]; - next if ($r1sam_flag == 4); - - # 77 = R1, unmapped - # 141 = R2, unmapped - - # 99 = R1, mapped, fwd strand --- frag on rev strand - # 101 = R1 unmapped, R2 mapped rev strand -- frag on rev strand - # 73 = R1, mapped, fwd strand --- frag on rev strand - # 147 = R2, mapped, rev strand -- frag on rev strand - # 153 = R2 mapped (R1 unmapped), rev strand -- frag on rev strand - # 133 = R2 unmapped, R1 mapped fwd strand -- frag on rev strand - - # 83 = R1, mapped, rev strand --- frag on fwd strand - # 69 = R1 unmapped, R2 mapped fwd strand -- frag on fwd strand - # 89 = R1 mapped rev strand, R2 unmapped -- frag on fwd strand - # 163 = R2, mapped, fwd strand -- frag on fwd strand - # 137 = R2 mapped (R1 unmapped), fwd strand -- frag on fwd strand - # 165 = R2 unmapped, R1 rev strand -- frag on fwd strand - - - my $frag_strand; -### This section is for only properly paired reads -# if ($r1sam_flag == 99) { -# $frag_strand = "-"; -# } elsif ($r1sam_flag == 83) { -# $frag_strand = "+"; -# } elsif ($r1sam_flag == 147) { -# $frag_strand = "-"; -# ($r1,$r2) = ($r2,$r1); -# @tmp_r1 = split(/\t/,$r1); -# @tmp_r2 = split(/\t/,$r2); -# } elsif ($r1sam_flag == 163) { -# ($r1,$r2) =($r2,$r1); -# @tmp_r1 = split(/\t/,$r1); -# @tmp_r2 = split(/\t/,$r2); -# $frag_strand = "+"; -# } else { -# print STDERR "R1 strand error $r1sam_flag\n"; -# next; -# } -### - if ($r1sam_flag == 16 || $r1sam_flag == 272) { - $frag_strand = "-"; - } elsif ($r1sam_flag eq "0" || $r1sam_flag == 256) { - $frag_strand = "+"; - } else { -# print STDERR "R1 strand error $r1sam_flag\n"; - next; - } - - - my @read_name = split(/\:/,$tmp_r1[0]); - - my $r1_cigar = $tmp_r1[5]; - my $r1_chr = $tmp_r1[2]; - my $r1_start = $tmp_r1[3]; - my $flags_r1 = join("\t",@tmp_r1[11..$#tmp_r1]); - - my $r1_mismatch; - if ($flags_r1 =~ /MD\:Z\:(\S+)\s/ || $flags_r1 =~ /MD\:Z\:(\S+?)$/) { - $r1_mismatch = $1; - } - - my $r1_phred = $tmp_r1[10]; - my $r1_seq = $tmp_r1[9]; - my $r1_mmscore = &get_alignment_score($r1_mismatch,$r1_cigar,$r1_phred,$r1_seq); - my $total_mmscore = $r1_mmscore; - -#check if read already exists from family mapping - 8/31/16 - this runs REGARDLESS of whether unique mapping overlaps a RepBase element or not - if (exists $read_hash{$r1name}{file1flag}) { - unless (exists $read_hash{$r1name}{rep_score}) { - print STDERR "err $r1name\n"; - } - if ($total_mmscore > $read_hash{$r1name}{rep_score} + 24) { - - # if unique mapping to genome is more than 2 mismatches per read = 2 * 2 * 6 alignment score better than to repeat element, throw out repeat element and use genome mapping -# print "deleting $r1name from rep element mapping - $total_mmscore $read_hash{$r1name}{rep_score}\n$r1\n$read_hash{$r1name}{R1}\n----\n"; - $deleted{$read_hash{$r1name}{ensttype}}++; - delete($read_hash{$r1name}) if (exists $read_hash{$r1name}); - - } else { - # repeat mapping exists and is better than unique genome mapping - throw out unique genome mapping - next; - } - } - - - - my @read_regions = &parse_cigar_string($r1_start,$r1_cigar,$r1_chr,$frag_strand); - -#this section is doing overlaps versus "peaks" = unique genomic repetitive elements = bed file regions from the repmasker run on hg19 - my %tmp_hash; - for my $region (@read_regions) { - my ($rchr,$rstr,$rpos) = split(/\:/,$region); - my ($rstart,$rstop) = split(/\-/,$rpos); - - my $verbose_flag = 0; - - my $rx = int($rstart / $genome_hashing_value); - my $ry = int($rstop / $genome_hashing_value); - for my $ri ($rx..$ry) { - - for my $peak (@{$peaks{$rchr}{$rstr}{$ri}}) { - my ($pchr,$ppos,$pstr,$ptype) = split(/\:/,$peak); - my ($pstart,$pstop) = split(/\-/,$ppos); - -# next if ($pstart > $rstop || $pstop < $rstart); -# if ($pstop == $rstart) { -# if ($pstart == $rstop) { -# if ($pstart == $rstart) { -# print "region $region $r1 peak $peak\n"; -# } - next if ($pstart >= $rstop || $pstop <= $rstart); - $tmp_hash{$peak} = 1; - - } - } - } - - my %temp_peak_read_counts; - my %converted_types; - for my $peak (keys %tmp_hash) { - my ($pchr,$ppos,$pstr,$ptype) = split(/\:/,$peak); - $temp_peak_read_counts{$ptype."_uniquegenomic"}++; - $converted_types{$convert_enst2type{$ptype}} = 1; - print STDERR "peak $peak $ptype\n" unless (exists $convert_enst2type{$ptype}); - } - - my @sorted_types = sort {$a cmp $b} keys %temp_peak_read_counts; - my $all_mapped_ensts = join("|",@sorted_types); - my @sorted_convertedtype = sort {$a cmp $b} keys %converted_types; - my $all_ensttypes = join("|",@sorted_convertedtype); - - -## 20191023 - in my re-parsing, I assign anything that's in a microRNA or microRNA proximal region to those classes regardless of whether they overlap a repbase element or not - copying that to this script for now - - - if (exists $converted_types{"miRNA"} || exists $converted_types{"miRNA-proximal"}) { - my $final_feature_type; - if (exists $converted_types{"miRNA"}) { - $final_feature_type = "miRNA"; - } elsif (exists $converted_types{"miRNA-proximal"}) { - $final_feature_type = "miRNA-proximal"; - } - - $all_ensttypes = "unique_".$final_feature_type; -#20171019 - changed this to sort to be consistent across perl versions -# my @sorted_keys = sort {$a cmp $b} keys %{$tmp_mirbase_hash{$final_feature_type}}; -# $all_mapped_ensts = "unique_".join("|",@sorted_keys); - } else { - - # now throw out reads that uniquely map but don't overlap a RepBase element - removed for now, can add back in later - # clarifying note: what I'm doing here is if it overlaps with >=1 repbase element, it gets counted as THAT; otherwise, it goes through and gets counted as CDS/proxintron/etc - unless (scalar(keys %temp_peak_read_counts) >= 1) { - - my $read1_start_position; - if ($frag_strand eq "+") { - $read1_start_position = $r1_start; - } elsif ($frag_strand eq "-") { - my $last_region = $read_regions[$#read_regions]; - my ($rchr,$rstr,$rpos) = split(/\:/,$last_region); - my ($rstart,$rstop) = split(/\-/,$rpos); - $read1_start_position = $rstop - 1; - } else { - print STDERR "error $frag_strand\n"; - } - - my $feature_flag = 0; - my %tmp_gencode_hash; - - my $rx = int($read1_start_position / $genome_hashing_value); - for my $gencode (@{$gencode_features{$r1_chr}{$frag_strand}{$rx}}) { - my ($gencode_enst,$gencode_type,$gencode_region) = split(/\|/,$gencode); - my ($gencode_start,$gencode_stop) = split(/\-/,$gencode_region); - - next if ($read1_start_position < $gencode_start); - next if ($read1_start_position >= $gencode_stop); - my $gencode_ensg = $enst2ensg{$gencode_enst}; - $tmp_gencode_hash{$gencode_type}{$gencode_ensg}="contained"; - $feature_flag = 1; - - } - - my $final_feature_type = "intergenic"; - - if ($feature_flag == 1) { - if (exists $tmp_gencode_hash{"CDS"}) { - $final_feature_type = "CDS"; - } elsif (exists $tmp_gencode_hash{"3utr"} || exists $tmp_gencode_hash{"5utr"}) { - if (exists $tmp_gencode_hash{"3utr"} && exists $tmp_gencode_hash{"5utr"}) { - $final_feature_type = "5utr_and_3utr"; - } elsif (exists $tmp_gencode_hash{"3utr"}) { - $final_feature_type = "3utr"; - } elsif (exists $tmp_gencode_hash{"5utr"}) { - $final_feature_type = "5utr"; - } else { - print STDERR "weird shouldn't hit this\n"; - } - } elsif (exists $tmp_gencode_hash{"proxintron"}) { - $final_feature_type = "proxintron"; - } elsif (exists $tmp_gencode_hash{"distintron"}) { - $final_feature_type = "distintron"; - } elsif (exists $tmp_gencode_hash{"noncoding_exon"}) { - $final_feature_type = "noncoding_exon"; - } elsif (exists $tmp_gencode_hash{"noncoding_proxintron"}) { - $final_feature_type = "noncoding_proxintron"; - } elsif (exists $tmp_gencode_hash{"noncoding_distintron"}) { - $final_feature_type = "noncoding_distintron"; - } elsif (exists $tmp_gencode_hash{"antisense_gencode"}) { - $final_feature_type = "antisense_gencode"; - } else { - print STDERR "weird - shouldn't hit this? elements of tmp_gencode_hash are ".join("|",keys %tmp_gencode_hash)."\n"; - } - } - -# $read_counts{$final_feature_type}++; - - - - $all_ensttypes = "unique_".$final_feature_type; - if (exists $tmp_gencode_hash{$final_feature_type}) { - -#20171019 - changed this to sort to be consistent across perl versions - my @sorted_keys = sort {$a cmp $b} keys %{$tmp_gencode_hash{$final_feature_type}}; - $all_mapped_ensts = "unique_".join("|",@sorted_keys); -# $all_mapped_ensts = "unique_".join("|",keys %{$tmp_gencode_hash{$final_feature_type}}); - } else { - $all_mapped_ensts = "unique_".$final_feature_type; - } - } - } -# next unless (scalar(keys %temp_peak_read_counts) >= 1); - - - -# print "keeping new mapping $all_mapped_ensts \t $all_ensttypes\n$r1\n"; - - - $read_hash{$r1name}{rep_flag} = 0; - $read_hash{$r1name}{rep_score} = $total_mmscore; - $read_hash{$r1name}{file2flag} = 1; - $read_hash{$r1name}{R1} = $r1; - $read_hash{$r1name}{mult_ensts} = $all_mapped_ensts; - $read_hash{$r1name}{ensttype} = $all_ensttypes; - - } - close(B); -} - - -sub read_rep_family_pe { - my $sam_file_familymapping = shift; -# print STDERR "opening $sam_file_familymapping\n"; - open(BAM,$sam_file_familymapping) || die "can't open $sam_file_familymapping\n"; -#open(BAM,"-|", "samtools view $rmDup_bam"); \ - - - while () { - my $r1 = $_; - chomp($r1); - if ($r1 =~ /^\@/) { - next; - } - - my $r2 = ; - chomp($r2); - - my @tmp_r1 = split(/\t/,$r1); - my @tmp_r2 = split(/\t/,$r2); - - my ($r1name,$r1bc) = split(/\s+/,$tmp_r1[0]); - my ($r2name,$r2bc) = split(/\s+/,$tmp_r2[0]); - - unless ($r1name eq $r2name) { - print STDERR "paired end mismatch error: $sam_file_familymapping r1 $tmp_r1[0] r2 $tmp_r2[0]\n"; - } - - my $r1sam_flag = $tmp_r1[1]; - my $r2sam_flag = $tmp_r2[1]; - unless ($r1sam_flag) { - print STDERR "error $r1 $r2\n"; - } - next if ($r1sam_flag == 77 || $r1sam_flag == 141); - - my $frag_strand; -### This section is for only properly paired reads - if ($r1sam_flag == 99 || $r1sam_flag == 355) { - $frag_strand = "-"; - } elsif ($r1sam_flag == 83 || $r1sam_flag == 339) { - $frag_strand = "+"; - } elsif ($r1sam_flag == 147 || $r1sam_flag == 403) { - ($r1,$r2) =($r2,$r1); - @tmp_r1 = split(/\t/,$r1); - @tmp_r2 = split(/\t/,$r2); - $frag_strand = "-"; - } elsif ($r1sam_flag == 163 || $r1sam_flag == 419) { - ($r1,$r2) =($r2,$r1); - @tmp_r1 = split(/\t/,$r1); - @tmp_r2 = split(/\t/,$r2); - $frag_strand = "+"; - } else { - next; - print STDERR "R1 strand error $r1sam_flag\n"; - } - - my $flags_r1 = join("\t",@tmp_r1[11..$#tmp_r1]); - my $flags_r2 = join("\t",@tmp_r2[11..$#tmp_r2]); - - my $r1_score; - if ($flags_r1 =~ /AS\:i\:(\S+)\s/ || $flags_r1 =~ /AS\:i\:(\S+)$/) { - $r1_score = $1; - } else { - print STDERR "couldn't find score for r1 $r1\n"; - } - my $r2_score; - if ($flags_r2 =~ /AS\:i\:(\S+)\s/ || $flags_r2 =~ /AS\:i\:(\S+)$/) { - $r2_score = $1; - } else { - print STDERR "couldn't find score for r2 $r2\n"; - } - my $total_score = $r1_score + $r2_score; - - my $all_mapped_ensts; - if ($flags_r1 =~ /ZZ\:Z\:(\S+?)\s/ || $flags_r1 =~ /ZZ\:Z\:(\S+?)$/) { - $all_mapped_ensts = $1; - } else { - print STDERR "didn't match? $flags_r1\n"; - } - - # if read has never been seen before, keep first mapping - my $all_info = $tmp_r1[2]; - my ($all_ensttypes,$all_primary_enst) = split(/\|\|/,$all_info); -# my $all_primary_enst = $tmp_r1[2]; - - my @ensttypes = split(/\|/,$all_ensttypes); - my @mapped_ensts = split(/\|/,$all_primary_enst); - $read_hash{$r1name}{rep_flag} = 1; - $read_hash{$r1name}{rep_score} = $total_score; - $read_hash{$r1name}{file1flag} = 1; - $read_hash{$r1name}{R1} = $r1; - $read_hash{$r1name}{R2} = $r2; - - my $sorted_mapped_ensts = &sort_priority($all_mapped_ensts); - $read_hash{$r1name}{mult_ensts} = $sorted_mapped_ensts; - $read_hash{$r1name}{ensttype} = $all_ensttypes; - } - close(BAM); - -} - - -sub read_rep_family_se { - - my $sam_file_familymapping = shift; -# print STDERR "opening $sam_file_familymapping\n"; - open(BAM,$sam_file_familymapping) || die "can't open $sam_file_familymapping\n"; -#open(BAM,"-|", "samtools view $rmDup_bam"); - - while () { - my $r1 = $_; - chomp($r1); - if ($r1 =~ /^\@/) { - next; - } - - my @tmp_r1 = split(/\t/,$r1); - my ($r1name,$r1bc) = split(/\s+/,$tmp_r1[0]); - my $r1sam_flag = $tmp_r1[1]; -# next if ($r1sam_flag == 77 || $r1sam_flag == 141); - next if ($r1sam_flag == 4); - - - my $frag_strand; - if ($r1sam_flag == 16 || $r1sam_flag == 272) { - $frag_strand = "-"; - } elsif ($r1sam_flag eq "0" || $r1sam_flag == 256) { - $frag_strand = "+"; - } else { - next; - print STDERR "R1 strand error $r1sam_flag\n"; - } - - -### This section is for only properly paired reads -# if ($r1sam_flag == 99 || $r1sam_flag == 355) { -# $frag_strand = "-"; -# } elsif ($r1sam_flag == 83 || $r1sam_flag == 339) { -# $frag_strand = "+"; -# } elsif ($r1sam_flag == 147 || $r1sam_flag == 403) { - # ($r1,$r2) =($r2,$r1); -# @tmp_r1 = split(/\t/,$r1); -# @tmp_r2 = split(/\t/,$r2); -# $frag_strand = "-"; -# } elsif ($r1sam_flag == 163 || $r1sam_flag == 419) { - # ($r1,$r2) =($r2,$r1); -# @tmp_r1 = split(/\t/,$r1); -# @tmp_r2 = split(/\t/,$r2); -# $frag_strand = "+"; -# } else { -# next; -# print STDERR "R1 strand error $r1sam_flag\n"; -# } -### - my $flags_r1 = join("\t",@tmp_r1[11..$#tmp_r1]); - - my $r1_score; - if ($flags_r1 =~ /AS\:i\:(\S+)\s/ || $flags_r1 =~ /AS\:i\:(\S+)$/) { - $r1_score = $1; - } else { - print STDERR "couldn't find score for r1 $r1\n"; - } - my $total_score = $r1_score; - - my $all_mapped_ensts; - if ($flags_r1 =~ /ZZ\:Z\:(\S+?)\s/ || $flags_r1 =~ /ZZ\:Z\:(\S+?)$/) { - $all_mapped_ensts = $1; - } else { - print STDERR "didn't match? $flags_r1\n"; - } - - # if read has never been seen before, keep first mapping - my $all_info = $tmp_r1[2]; - my ($all_ensttypes,$all_primary_enst) = split(/\|\|/,$all_info); -# my $all_primary_enst = $tmp_r1[2]; - my @ensttypes = split(/\|/,$all_ensttypes); - my @mapped_ensts = split(/\|/,$all_primary_enst); - $read_hash{$r1name}{rep_flag} = 1; - $read_hash{$r1name}{rep_score} = $total_score; - $read_hash{$r1name}{file1flag} = 1; - $read_hash{$r1name}{R1} = $r1; - - my $sorted_mapped_ensts = &sort_priority($all_mapped_ensts); - $read_hash{$r1name}{mult_ensts} = $sorted_mapped_ensts; - -# $read_hash{$r1name}{mult_ensts} = $all_mapped_ensts; - - $read_hash{$r1name}{ensttype} = $all_ensttypes; - } - close(BAM); - -} - - - -sub read_peakfi { - my $fi = shift; - my $duplicate = 0; - my $lc=0; -# print STDERR "reading peak file $fi\n"; - open(F,$fi) || die "no peak file $fi\n"; - while () { - my $line = $_; - chomp($line); - my @tmp = split(/\t/,$line); - my $chr = shift(@tmp); - my $start = shift(@tmp); - my $stop = shift(@tmp); - my $gene = shift(@tmp); - my $pval = shift(@tmp); - my $strand = shift(@tmp); - -# $gene =~ s/\_$//g; - $gene = uc($gene); - - next if ($chr eq "genoName"); - - my $x = int($start / $genome_hashing_value); - my $y = int($stop / $genome_hashing_value); - - my $peak = $chr.":".$start."-".$stop.":".$strand.":".$gene; -# if (exists $peak_read_counts{$peak}{allpeaks}) { -# $duplicate++; -# } else { -# $lc++; -# } -# $peak_read_counts{$gene}{allpeaks} = 1; -# $peak_read_counts{$gene."_antisense"}{allpeaks} = 1; - - - for my $i ($x..$y) { - push @{$peaks{$chr}{$strand}{$i}},$chr.":".$start."-".$stop.":".$strand.":".$gene; - push @{$peaks{$chr}{$revstrand{$strand}}{$i}},$chr.":".$start."-".$stop.":".$strand.":antisense_".$gene; - } - } - close(F); -# print STDERR "duplciate $duplicate\nlc $lc\n"; -} - -sub sort_priority { - my $enstlist_unsorted = shift; - my @enstlist_split = split(/\|/,$enstlist_unsorted); - - for my $enst (@enstlist_split) { - print "$enst\t".$convert_enst2priorityN{lc($enst)}."\n" unless (exists $convert_enst2priorityN{lc($enst)}); - } - my @sorted_enstlist = sort {$convert_enst2priorityN{lc($a)} <=> $convert_enst2priorityN{lc($b)}} @enstlist_split; - my $joined_sorted_enstlist = join("\|",@sorted_enstlist); - return($joined_sorted_enstlist); -} - -sub read_in_filelists { - my $fi = shift; - my $priority_n = 0; - my %corrected_ensts; - open(F,$fi); - for my $line () { - chomp($line); - my ($allenst,$allensg,$gid,$type_label,$typefile) = split(/\t/,$line); -# my ($allensg,$gid,$allenst) = split(/\t/,$line); - unless ($allenst) { - print STDERR "error missing enst $line $fi\n"; - } - $allenst = uc($allenst); - my @ensts = split(/\|/,$allenst); - $gid =~ s/\?$//; - $gid =~ s/\_$//; - $type_label =~ s/\?$//; - for my $enst (@ensts) { - if ($enst =~ /\_$/) { - $corrected_ensts{$enst}++; -# $enst =~ s/\_$//; - } - $enst2gene{$enst} = $gid; - $enst2gene{"antisense_".$enst} = "antisense_".$gid; -# $convert_enst2type{$enst} = $type_label.":".$priority_n; -# $convert_enst2type{$enst."_antisense"} = $type_label."_antisense:".$priority_n; - $convert_enst2type{$enst} = $type_label; - $convert_enst2type{"antisense_".$enst} = "antisense_".$type_label; - $convert_enst2priority{lc($enst)} = $type_label."|".$priority_n; - $convert_enst2priorityN{lc($enst)} = $priority_n; - $convert_enst2priorityN{lc($enst."_DOUBLEMAP")} = $priority_n+0.25; - $convert_enst2priorityN{"antisense_".lc($enst)} = $priority_n+0.5; - $convert_enst2priorityN{"antisense_".lc($enst."_DOUBLEMAP")} = $priority_n+0.75; - $convert_enst2priority{"antisense_".lc($enst)} = "antisense_".$type_label."|".$priority_n; - - $priority_n++; - } - } - close(F); - # for my $k (keys %corrected_ensts) { -# print "$k\t$corrected_ensts{$k}\n"; -# } -} - -sub min { - my $x = shift; - my $y = shift; - if ($x < $y) { - return($x); - } else { - return($y); - } -} - - - -sub get_alignment_score { - my $md_string = shift; - my $cigar_string = shift; - my $quality_string = shift; - my $read_seq = shift; - my (@insertions) = &parse_cigar_string_foralignmentscore($cigar_string,$quality_string); - my $gap_score = shift(@insertions); - my $mm_score = &parse_mismatch_string_foralignmentscore($md_string,$quality_string,\@insertions,$read_seq); - my $total_score = 0 - $mm_score - $gap_score; - return($total_score); -} - - - - -sub parse_cigar_string { - my $region_start_pos = shift; - my $flags = shift; - my $chr = shift; - my $strand = shift; - - my $current_pos = $region_start_pos; - my @regions; - - while ($flags =~ /(\d+)([A-Z])/g) { - - if ($2 eq "N") { - #read has intron of N bases at location - -# 1 based, closed ended to 0 based, right side open ended fix -# push @regions,$chr.":".$strand.":".$region_start_pos."-".($current_pos-1); - push @regions,$chr.":".$strand.":".($region_start_pos-1)."-".($current_pos-1); - - $current_pos += $1; - $region_start_pos = $current_pos; - } elsif ($2 eq "M") { - #read and genome match - $current_pos += $1; - } elsif ($2 eq "S") { - #beginning of read is soft-clipped; mapped pos is actually start pos of mapping not start of read - } elsif ($2 eq "I") { - #read has insertion relative to genome; doesn't change genome position - } elsif ($2 eq "D") { -# push @read_regions,$chr.":".$current_pos."-".($current_pos+=$1); - $current_pos += $1; - #read has deletion relative to genome; genome position has to increase - } else { - print STDERR "flag $1 $2 $flags\n"; - - } - } - -# 1 based, closed ended to 0 based, right side open ended fix -# $region_start_pos is 1-based, closed ended -> ($region_start_pos-1) is 0-based, closed ended -# ($current_pos-1) is 1-based, closed ended -> ($current_pos-1-1) is 0-based, closed ended -> ($current_pos-1-1+1) is 0-based, open ended - -# push @regions,$chr.":".$strand.":".$region_start_pos."-".($current_pos-1); - push @regions,$chr.":".$strand.":".($region_start_pos-1)."-".($current_pos-1); - - return(@regions); -} - - - - - -sub parse_mismatch_string_foralignmentscore { - my $flags = shift; - my $phred_scores = shift; - my $insertionsref = shift; - my $read_seq = shift; - - my @insertions = @$insertionsref; - my ($next_insertion_pos,$next_insertion_len); - my $insertion_flag = 0; - if (@insertions > 0) { - if ($insertions[0] eq "none") { - $insertion_flag = 0; - } else { - $insertion_flag = 1; - my $next_insertion = shift(@insertions); - ($next_insertion_pos,$next_insertion_len) = split(/\|/,$next_insertion); - } - } else { - $insertion_flag = 0; - } - - if ($flags =~ /^\d+$/) { - return(0); - } - my @phred_scoress = split(//,$phred_scores); - my @quality_scores; - for (my $i=0;$i<@phred_scoress;$i++) { - $quality_scores[$i] = ord($phred_scoress[$i]) - 33; - print STDERR "error quality score out of range $quality_scores[$i] $phred_scoress[$i]\n" if ($quality_scores[$i] < 0 || $quality_scores[$i] > 43); -# $quality_scores[$i] = $convert_phred{$phred_scoress[$i]}; - } - - my $current_pos = 0; - - my $mm_penalty = 0; - my @regions; - my $mn = 2; - my $mx = 6; - - my $current_read_pos = 0; - while ($flags) { - if ($flags =~ /^(\d+)/) { - $current_read_pos += $1; - if ($insertion_flag == 1) { - while ($insertion_flag == 1 && $current_read_pos >= $next_insertion_pos) { - $current_read_pos += $next_insertion_len; - if (@insertions > 0) { - $insertion_flag = 1; - my $next_insertion = shift(@insertions); - ($next_insertion_pos,$next_insertion_len) = split(/\|/,$next_insertion); - } else { - $insertion_flag = 0; - } - } - } - - $flags = substr($flags,length($1)); - } elsif ($flags =~ /^([A-Z]+)/) { - my $len = length($1); - - for my $j (1..$len) { - my $base_score = $quality_scores[$current_read_pos+$j-1]; - unless ($base_score) { - print STDERR "$j $j cur $current_read_pos read $read_seq $phred_scores\n"; - } - my $base = substr($1,$j-1,1); - my $base_seq = substr($read_seq,$current_read_pos+$j-1,1); - my $base_mm_penalty; - if ($base eq "N" || $base_seq eq "N") { - $base_mm_penalty = 1; - } elsif (($base eq "A" || $base eq "C" || $base eq "G" || $base eq "T") && ($base_seq eq "A" || $base_seq eq "C" || $base_seq eq "G" || $base_seq eq "T")) { - $base_mm_penalty = $mn + floor( ($mx-$mn) * (&min($base_score,40)/40) ); - } else { - print STDERR "unexpected base - $base\n"; - } - $mm_penalty += $base_mm_penalty; - } - $current_read_pos += $len; - $flags = substr($flags,length($1)); - } elsif ($flags =~ /^\^([A-Z]+)\d/) { - my $len = length($1); - #insertions aren't penalized here; penalized in gap open/close - $flags = substr($flags,$len+1); - } else { - print STDERR "this is a flag I'm not expecting $flags\n"; - } - } - return($mm_penalty); -} - -sub parse_cigar_string_foralignmentscore { - my $flags = shift; - my $phred_scores = shift; - - if ($flags =~ /^\d+M$/) { - return(0); - } - my @phred_scoress = split(//,$phred_scores); - my @quality_scores; - for (my $i=0;$i<@phred_scoress;$i++) { - $quality_scores[$i] = ord($phred_scoress[$i]) - 33; - print STDERR "error quality score out of range $quality_scores[$i] $phred_scoress[$i]\n" if ($quality_scores[$i] < 0 || $quality_scores[$i] > 43); -# $quality_scores[$i] = $convert_phred{$phred_scoress[$i]}; - } - - my $gap_open = 5; - my $gap_extend = 3; - my $current_pos = 0; - my $mm_penalty = 0; - my @insertions; - my $current_read_pos = 0; - - while ($flags =~ /(\d+)([A-Z])/g) { - if ($2 eq "N") { - # intron - not sure what to do with this for now? I think skip - - $current_read_pos += $1; - $current_pos += $1; - } elsif ($2 eq "M") { - $current_read_pos += $1; - $current_pos += $1; - } elsif ($2 eq "S") { - $current_read_pos += $1; - $mm_penalty += $gap_open + $1 * $gap_extend; - } elsif ($2 eq "I") { - push @insertions,$current_read_pos."|".$1; - $current_read_pos += $1; - $mm_penalty += $gap_open + $1 * $gap_extend; - } elsif ($2 eq "D") { - $current_pos += $1; - $mm_penalty += $gap_open + $1 * $gap_extend; - } else { - print STDERR "flag $1 $2 $flags\n"; - - } - } - return($mm_penalty,@insertions); - -} - - -sub read_gencode { - ## eric note: this has been tested for off-by-1 issues with ucsc brower table output! \ - - my $fi = shift; -# my $fi = "/projects/ps-yeolab/genomes/hg19/gencode_v19/gencodev19_comprehensive"; - print STDERR "reading in $fi\n"; - open(F,$fi) || die "couldn't open gencode $fi\n"; - while () { - chomp($_); - my @tmp = split(/\t/,$_); - my $enst = $tmp[1]; - next if ($enst eq "name"); - my $chr = $tmp[2]; - my $str = $tmp[3]; - my $txstart = $tmp[4]; - my $txstop = $tmp[5]; - my $cdsstart = $tmp[6]; - my $cdsstop = $tmp[7]; - - my @starts = split(/\,/,$tmp[9]); - my @stops = split(/\,/,$tmp[10]); - - my @tmp_features; - - my $transcript_type = $enst2type{$enst}; - unless ($transcript_type) { - print STDERR "error transcript_type $transcript_type $enst\n"; - } - if ($transcript_type eq "protein_coding") { - - for (my $i=0;$i<@starts;$i++) { - if ($str eq "+") { - if ($stops[$i] < $cdsstart) { - # exon is all 5' utr \ - - push @tmp_features,$enst."|5utr|".($starts[$i])."-".$stops[$i]; - } elsif ($starts[$i] > $cdsstop) { - #exon is all 3' utr \ - - push @tmp_features,$enst."|3utr|".($starts[$i])."-".$stops[$i]; - } elsif ($starts[$i] > $cdsstart && $stops[$i] < $cdsstop) { - #exon is all coding \ - - push @tmp_features,$enst."|CDS|".($starts[$i])."-".$stops[$i]; - } else { - my $cdsregion_start = $starts[$i]; - my $cdsregion_stop = $stops[$i]; - - if ($starts[$i] <= $cdsstart && $cdsstart <= $stops[$i]) { - #cdsstart is in exon \ - - my $five_region = ($starts[$i])."-".$cdsstart; - push @tmp_features,$enst."|5utr|".$five_region; - $cdsregion_start = $cdsstart; - } - - if ($starts[$i] <= $cdsstop && $cdsstop <= $stops[$i]) { - #cdsstop is in exon \ - - my $three_region = ($cdsstop)."-".$stops[$i]; - push @tmp_features,$enst."|3utr|".$three_region; - $cdsregion_stop = $cdsstop; - } - my $cds_region = ($cdsregion_start)."-".$cdsregion_stop; - push @tmp_features,$enst."|CDS|".$cds_region; - } - } elsif ($str eq "-") { - if ($stops[$i] < $cdsstart) { - # exon is all 5' utr \ - - push @tmp_features,$enst."|3utr|".($starts[$i])."-".$stops[$i]; - } elsif ($starts[$i] > $cdsstop) { - #exon is all 3' utr \ - - push @tmp_features,$enst."|5utr|".($starts[$i])."-".$stops[$i]; - } elsif ($starts[$i] > $cdsstart &&$stops[$i] < $cdsstop) { - #exon is all coding \ - - push @tmp_features,$enst."|CDS|".($starts[$i])."-".$stops[$i]; - } else { - my $cdsregion_start = $starts[$i]; - my $cdsregion_stop = $stops[$i]; - - if ($starts[$i] <= $cdsstart && $cdsstart <= $stops[$i]) { - #cdsstart is in exon \ - - my $three_region = ($starts[$i])."-".$cdsstart; - push @tmp_features,$enst."|3utr|".$three_region; - $cdsregion_start = $cdsstart; - } - - if ($starts[$i] <= $cdsstop && $cdsstop <= $stops[$i]) { - #cdsstop is in exon \ - - my $five_region = ($cdsstop)."-".$stops[$i]; - push @tmp_features,$enst."|5utr|".$five_region; - $cdsregion_stop = $cdsstop; - } - - my $cds_region = ($cdsregion_start)."-".$cdsregion_stop; - push @tmp_features,$enst."|CDS|".$cds_region; - } - } - } - for (my $i=0;$i 2 * 500) { - push @tmp_features,$enst."|proxintron|".($stops[$i])."-".($stops[$i]+500); - push @tmp_features,$enst."|distintron|".($stops[$i]+500)."-".($starts[$i+1]-500); - push @tmp_features,$enst."|proxintron|".($starts[$i+1]-500)."-".$starts[$i+1]; - } else { - my $midpoint = int(($starts[$i+1]+$stops[$i])/2); - push @tmp_features,$enst."|proxintron|".($stops[$i])."-".($midpoint); - push @tmp_features,$enst."|proxintron|".($midpoint)."-".$starts[$i+1]; - } -## push @tmp_features,$enst."|intron|".($stops[$i])."-".$starts[$i+1]; - } - } else { - - for (my $i=0;$i<@starts;$i++) { - push @tmp_features,$enst."|noncoding_exon|".($starts[$i])."-".$stops[$i]; - } - for (my $i=0;$i 2 * 500) { - push @tmp_features,$enst."|noncoding_proxintron|".($stops[$i])."-".($stops[$i]+500); - push @tmp_features,$enst."|noncoding_distintron|".($stops[$i]+500)."-".($starts[$i+1]-500); - push @tmp_features,$enst."|noncoding_proxintron|".($starts[$i+1]-500)."-".$starts[$i+1]; - } else { - my $midpoint = int(($starts[$i+1]+$stops[$i])/2); - push @tmp_features,$enst."|noncoding_proxintron|".($stops[$i])."-".($midpoint); - push @tmp_features,$enst."|noncoding_proxintron|".($midpoint)."-".$starts[$i+1]; - } - - -# push @tmp_features,$enst."|noncoding_intron|".($stops[$i])."-".$starts[$i+1]; - } - } - - for my $feature (@tmp_features) { - my ($enst,$type,$region) = split(/\|/,$feature); - my ($reg_start,$reg_stop) = split(/\-/,$region); - my $x = int($reg_start/$genome_hashing_value); - my $y = int($reg_stop /$genome_hashing_value); - - for my $j ($x..$y) { - push @{$gencode_features{$chr}{$str}{$j}},$feature; - push @{$gencode_features{$chr}{$revstrand{$str}}{$j}},$enst."|antisense_gencode|".$reg_start."-".$reg_stop; - } - } - } - close(F); - -} - - - -sub read_gencode_gtf { - - my $file = shift; -# my $file = "/projects/ps-yeolab/genomes/hg19/gencode_v19/gencode.v19.chr_patch_hapl_scaff.annotation.gtf"; - print STDERR "Reading in $file\n"; - open(F,$file) || die "couldn't open gencode $file\n"; - for my $line () { - chomp($line); - next if ($line =~ /^\#/); - my @tmp = split(/\t/,$line); - next unless ($tmp[2] eq "transcript"); - - my $stuff = $tmp[8]; - my @stufff = split(/\;/,$stuff); - my ($ensg_id,$gene_type,$gene_name,$enst_id,$transcript_type); - - for my $s (@stufff) { - $s =~ s/^\s//g; - $s =~ s/\s$//g; - - if ($s =~ /gene_id \"(.+?)\"/) { - if ($ensg_id) { - print STDERR "two ensg ids? $line\n"; - } - $ensg_id = $1; - } - if ($s =~ /transcript_id \"(.+?)\"/) { - if ($enst_id) { - print STDERR "two enst ids? $line\n"; - } - $enst_id = $1; - } - if ($s =~ /gene_type \"(.+?)\"/) { - if ($gene_type) { - print STDERR "two gene types $line\n"; - } - $gene_type = $1; - - } - - if ($s =~ /transcript_type \"(.+?)\"/) { - $transcript_type = $1; - } - if ($s =~ /gene_name \"(.+?)\"/) { - $gene_name = $1; - } - } - - if (exists $enst2ensg{$enst_id} && $ensg_id ne $enst2ensg{$enst_id}) { - print STDERR "error two ensgs for enst $enst_id ensgs $ensg_id $enst2ensg{$enst_id}\n$line\n"; - } - $enst2ensg{$enst_id} = $ensg_id; - $ensg2name{$ensg_id}{$gene_name}=1; - $ensg2type{$ensg_id}{$gene_type}=1; - $enst2type{$enst_id} = $transcript_type; - } - close(F); -} - - -sub get_type_flag { - my $ref = shift; - my %feature_hash = %$ref; - my $feature_type = shift; - - my $feature_type_final = "contained"; - for my $ensg (keys %{$feature_hash{$feature_type}}) { - $feature_type_final = "partial" unless ($feature_hash{$feature_type}{$ensg} eq "contained"); - } - return($feature_type_final); -} - - - - -sub reorder_transcripts_by_priority { - my $enst_ref = shift; - my @enst_list = @$enst_ref; - - my %reptype_hash; - for my $enst_id_orig (@enst_list) { - my $enst_id = $enst_id_orig; - $enst_id =~ s/\_DOUBLEMAP$//; - $enst_id =~ s/\_spliced$//; -# $enst_id =~ s/\_$//; - $enst_id = lc($enst_id); - unless (exists $convert_enst2priority{$enst_id} && $convert_enst2priority{$enst_id}) { - print STDERR "err $enst_id $enst_id_orig$convert_enst2priority{$enst_id}\n"; - } - my ($rep_type,$priority_n) = split(/\|/,$convert_enst2priority{$enst_id}); - $reptype_hash{$rep_type}{$enst_id_orig} = $priority_n; - } - - my @sorted_reptypes = sort {$a cmp $b} keys %reptype_hash; - my @sorted_ensts_bytype; - my @sorted_topenst_bytype; - for my $sorted_reptype (@sorted_reptypes) { - my @sorted_ensts_by_priority = sort {$reptype_hash{$sorted_reptype}{$a} <=> $reptype_hash{$sorted_reptype}{$b}} keys %{$reptype_hash{$sorted_reptype}}; - my $reptype_joined = join(";;",@sorted_ensts_by_priority); - push @sorted_ensts_bytype,$reptype_joined; - push @sorted_topenst_bytype,$sorted_ensts_by_priority[0]; - } - my $output_mult_ensts = join("|",@sorted_ensts_bytype); - my $output_top_ensts = join("|",@sorted_topenst_bytype); - my $output_types = join("|",@sorted_reptypes); - - my $primary_output = $output_types."||".$output_top_ensts; - my $full_output = $output_types."||".$output_mult_ensts; - return($primary_output,$full_output); -} - - diff --git a/bin/perl/merge_multiple_parsed_files.simplified_20191022.pl b/bin/perl/merge_multiple_parsed_files.simplified_20191022.pl deleted file mode 100755 index adcf0a2..0000000 --- a/bin/perl/merge_multiple_parsed_files.simplified_20191022.pl +++ /dev/null @@ -1,116 +0,0 @@ -#!/usr/bin/env perl - -use warnings; -use strict; -use Fcntl qw(:flock); - -my $output_fi = shift(@ARGV); -open(OUT,">$output_fi") || die "can't open output file $output_fi\n"; - -my @files = @ARGV; -my @split_fi1 = split(/\//,$files[0]); -my $short_fi1 = pop(@split_fi1); -my $working_dir = join("/",@split_fi1); -my $failed_jobs_list = $working_dir."/"."failed_jobs_list.txt"; - -print("SHORTFI:".$short_fi1); -print("WORKDIR:".$working_dir); - -my $lock_fi = $failed_jobs_list.".lck"; -open(S,">$lock_fi") || sleep 5; -flock(S, LOCK_EX) || sleep 5; -open(FH,">>$failed_jobs_list") or die "this is weird - shouldn't be hit $failed_jobs_list\n"; - -my %read_sums; -for my $file (@files) { -# my $donefi = $file.".done"; -# my $done_seq = `cat $donefi` if (-e $donefi); -# chomp($done_seq); -# if (-e $donefi) { -# print "donefi exists $donefi\n"; -# print "doneseq $done_seq X\n"; -# } - -# unless (-e $donefi && $done_seq eq "jobs done") { -# next; -# print FH "Failed job $output_fi on $file\n"; -# exit; -# } - &parse_file($file); -} -close(FH); -close(S); - -print OUT "#READINFO\tAllReads\t".$read_sums{all}."\n"; -print OUT "#READINFO\tUsableReads\t".$read_sums{usable}."\t".($read_sums{usable}/$read_sums{all})."\n"; -print OUT "#READINFO\tGenomicReads\t".$read_sums{genomic}."\t".($read_sums{genomic}/$read_sums{usable})."\n"; -print OUT "#READINFO\tRepFamilyReads\t".$read_sums{repfamily}."\t".($read_sums{repfamily}/$read_sums{usable})."\n"; - -my @sorted_total = sort {$read_sums{total}{$b} <=> $read_sums{total}{$a}} keys %{$read_sums{total}}; -for my $element (@sorted_total) { -#for my $element (keys %{$read_sums{total}}) { - print OUT "TOTAL\t$element\t".$read_sums{total}{$element}."\t".($read_sums{total}{$element}/$read_sums{usable})."\n"; -} - -my @sorted_element = sort {$read_sums{element}{$b}{readnum} <=> $read_sums{element}{$a}{readnum}} keys %{$read_sums{element}}; -for my $element (@sorted_element) { -#for my $element (keys %{$read_sums{element}}) { - print OUT "ELEMENT\t".$read_sums{element}{$element}{ensg_primary}."\t".$read_sums{element}{$element}{readnum}."\t".($read_sums{element}{$element}{readnum}/$read_sums{usable})."\t".$element."\t".$read_sums{element}{$element}{ensg_all}."\n"; -} - - - - -sub parse_file { - my $parsed_fi = shift; - open(P,$parsed_fi) || die "can't open $parsed_fi\n"; - while (

) { - chomp($_); - my @tmp = split(/\t/,$_); - if ($tmp[0] eq "#READINFO") { -# print STDERR "tmp0 $tmp[0] $_\n"; - my $line = $_; - if ($_ =~ /All\sreads\:\t(\d+)\tPCR\sduplicates\sremoved\:\t(\d+)\tUsable\sRemaining\:\t(\d+)\tUsable\sfrom\sgenomic\smapping\:\t(\d+)\tUsable\sfrom\sfamily\smapping\:\t(\d+)$/) { - - my $all_reads = $1; - my $pcr_duplicates = $2; - my $usable_reads = $3; - my $genomic_reads = $4; - my $repfamily_reads = $5; - - $read_sums{all} += $all_reads; -# $read_sums{usable} += $usable_reads; -# $read_sums{genomic} += $genomic_reads; -# $read_sums{repfamily} += $repfamily_reads; - -# print STDERR "all $all_reads rep $repfamily_reads genomic $genomic_reads usable $usable_reads\n"; - } elsif ($tmp[1] eq "UsableReads") { - $read_sums{usable} += $tmp[2]; - } elsif ($tmp[1] eq "GenomicReads") { - $read_sums{genomic} += $tmp[2]; - } elsif ($tmp[1] eq "RepFamilyReads") { - $read_sums{repfamily} += $tmp[2]; - } else { - print STDERR "couldn't parse readinfo line $_\n"; - } - } elsif ($tmp[0] eq "#READINFO2") { - } elsif ($tmp[0] eq "TOTAL") { - my ($total,$element,$readnum,$rpm) = @tmp; - - $read_sums{total}{$element} += $readnum; - } else { - my ($ele_delete,$ensg_primary,$readnum,$rpm,$enst_all,$ensg_all) = @tmp; - -# print "ELEMENT ".join("X x",@tmp)."\n"; - print STDERR "error - ensg_all mismatch $ensg_all $read_sums{element}{$enst_all}{ensg_all}\n" if (exists $read_sums{element}{$enst_all}{ensg_all} && ($read_sums{element}{$enst_all}{ensg_all} ne $ensg_all)); - print STDERR "error - ensg_primary mismatch $ensg_primary $read_sums{element}{$enst_all}{ensg_primary}\n" if (exists $read_sums{element}{$enst_all}{ensg_primary} && ($read_sums{element}{$enst_all}{ensg_primary} ne $ensg_primary)); - - $read_sums{element}{$enst_all}{ensg_all} = $ensg_all; - $read_sums{element}{$enst_all}{ensg_primary} = $ensg_primary; - $read_sums{element}{$enst_all}{readnum} += $readnum; - - - } - } - close(P); -} diff --git a/bin/perl/parse_bowtie2_output_realtime_includemultifamily_PE.pl b/bin/perl/parse_bowtie2_output_realtime_includemultifamily_PE.pl deleted file mode 100755 index 0017f3b..0000000 --- a/bin/perl/parse_bowtie2_output_realtime_includemultifamily_PE.pl +++ /dev/null @@ -1,520 +0,0 @@ -#!/usr/bin/env perl - -use warnings; -use strict; - -my %enst2gene; -my %convert_enst2type; -my %multimapping_hash; - -my %rRNA_extra_hash; -my %rRNA_extra_hash_rev; -$rRNA_extra_hash{"RNA28S"} = "RNA45S"; -$rRNA_extra_hash{"RNA18S"} = "RNA45S"; -$rRNA_extra_hash{"RNA5-8S"} = "RNA45S"; -$rRNA_extra_hash{"antisense_RNA28S"} = "antisense_RNA45S"; -$rRNA_extra_hash{"antisense_RNA18S"} = "antisense_RNA45S"; -$rRNA_extra_hash{"antisense_RNA5-8S"} = "antisense_RNA45S"; - -$rRNA_extra_hash_rev{"RNA45S"}{"RNA28S"} = 1; -$rRNA_extra_hash_rev{"RNA45S"}{"RNA18S"} = 1; -$rRNA_extra_hash_rev{"RNA45S"}{"RNA5-8S"} = 1; -$rRNA_extra_hash_rev{"antisense_RNA45S"}{"antisense_RNA28S"} = 1; -$rRNA_extra_hash_rev{"antisense_RNA45S"}{"antisense_RNA18S"} = 1; -$rRNA_extra_hash_rev{"antisense_RNA45S"}{"antisense_RNA5-8S"} = 1; - - - - -#my $working_directory = "/home/elvannostrand/data/clip/CLIPseq_analysis/ENCODE_v9_20151209/PolIII_mapping/"; -# my $working_directory = $ARGV[3]; -# my $species = $ARGV[5]; -my $filelist_file = $ARGV[4]; -#hg19 -# my $filelist_file = "/home/elvannostrand/data/clip/CLIPseq_analysis/RNA_type_analysis/MASTER_filelist.wrepbaseandtRNA.enst2id.fixed.UpdatedSimpleRepeat.20190424"; -# if ($species eq "hg38") { -# $filelist_file = "/home/elvannostrand/data/clip/CLIPseq_analysis/RNA_type_analysis/hg38/MASTER_FILELIST.20201203.wrepbaseandtRNA.enst2id.fixed.UpdatedSimpleRepeat.wmiRs.list"; -# $filelist_file = "/home/elvannostrand/data/clip/CLIPseq_analysis/RNA_type_analysis/hg38/MASTER_FILELIST.20201203.wrepbaseandtRNA.enst2id.fixed.UpdatedSimpleRepeat"; -# } -&read_in_filelists($filelist_file); - -my $print_batch = 10000; -my $read_counter = 0; - -#&read_in_filelists("/home/elvannostrand/data/clip/CLIPseq_analysis/RNA_type_analysis/genelists.RNA5S","RNA5S"); -#my $fastq_file1 = "/home/gpratt/projects/encode/analysis/encode_v9/484_CLIP_S6_L001_R1_001.B06_484_01_POLR2G.adapterTrim.round2.fastq.gz"; -#my $fastq_file2 = "/home/gpratt/projects/encode/analysis/encode_v9/484_CLIP_S6_L001_R2_001.B06_484_01_POLR2G.adapterTrim.round2.fastq.gz"; - -#my $fastq_file1 = "test_r1.fastq"; -#my $fastq_file2 = "test_r2.fastq"; -#my $fastq_file1 = "484_CLIP_S6_L001_R1_001.B06_484_01_POLR2G.adapterTrim.round2.fastq"; -#my $fastq_file2 = "484_CLIP_S6_L001_R2_001.B06_484_01_POLR2G.adapterTrim.round2.fastq"; -#my $fastq_file1 = "484_INPUT_S7_L001_R1_001.unassigned.adapterTrim.round2.fastq"; -#my $fastq_file2 = "484_INPUT_S7_L001_R2_001.unassigned.adapterTrim.round2.fastq"; - -my $fastq_file1 = $ARGV[0]; -my $fastq_file2 = $ARGV[1]; - -my $bowtie_db = $ARGV[2]; -#my $bowtie_db = "/home/elvannostrand/data/clip/CLIPseq_analysis/RNA_type_analysis/filelist_POLIII.combined.Nflank"; -my @bowtie_db_split = split(/\//,$bowtie_db); -my $bowtie_db_short = $bowtie_db_split[$#bowtie_db_split]; - -my @fastq_fi1_split = split(/\//,$fastq_file1); -my $fastq_fi1_short = $fastq_fi1_split[$#fastq_fi1_split]; - -my $output = $ARGV[3]; -#my $output = $working_directory.$fastq_fi1_short.".mapped_vs_".$bowtie_db_short.".sam"; -my $bowtie_out = $output.".bowtieout"; -open(SAMOUT,">$output"); -my $multimapping_out = $output.".multimapping_deleted"; -open(MULTIMAP,">$multimapping_out"); -my $done_file = $output.".done"; - -# -my $command = "stdbuf -oL bowtie2 -q --sensitive -a -p 1 --no-mixed --reorder -x $bowtie_db -1 $fastq_file1 -2 $fastq_file2 2> $bowtie_out"; -#my $command = "/projects/ps-yeolab/software/bowtie-1.1.1/./bowtie -q -1 $fastq_file1 -2 $fastq_file2 $bowtie_db -a -v 2 --best --strata -S 2> $bowtie_out"; -print STDERR "command $command\n"; -#my $pid = open(BOWTIE, "-|", "unbuffer /projects/ps-yeolab/software/bowtie-1.1.1/./bowtie -q -1 $fastq_file1 -2 $fastq_file2 $genelistfi -a -v 2 --best --strata -S 2> $bowtie_out"); -my $pid = open(BOWTIE, "-|", "stdbuf -oL bowtie2 -q --sensitive -a -p 3 --no-mixed --reorder -x $bowtie_db -1 $fastq_file1 -2 $fastq_file2 2> $bowtie_out"); -#my $pid = open(BOWTIE, "-|", "stdbuf -oL /projects/ps-yeolab/software/bowtie-1.1.1/./bowtie -q -1 $fastq_file1 -2 $fastq_file2 $bowtie_db -a -v 2 --best --strata -S 2> $bowtie_out"); -#print "PID $pid\n"; -my %fragment_hash; -my $duplicate_count=0; -my $unique_count=0; -my $all_count=0; -my $prev_r1name = ""; - -my %read_hash; -if ($pid) { - - while () { - my $r1 = $_; - chomp($r1); - if ($r1 =~ /^\@/) { - print SAMOUT "$r1\n"; - next; - } - - my $r2 = ; - chomp($r2); - - my @tmp_r1 = split(/\t/,$r1); - my @tmp_r2 = split(/\t/,$r2); - - my ($r1name,$r1bc) = split(/\s+/,$tmp_r1[0]); - my ($r2name,$r2bc) = split(/\s+/,$tmp_r2[0]); -# $all_count++; -# if ($all_count % 1000 == 0) { -# print STDERR "$all_count\n"; -# } - unless ($r1name eq $r2name) { - print STDERR "paired end mismatch error: r1 $tmp_r1[0] r2 $tmp_r2[0]\n"; - } - - my $debug_flag = 0; -# if ($r1name eq "GCCAGCCTTA:K00180:175:H7W2CBBXX:3:2106:25712:12075") { -# print "READ1\t$r1\n"; -# $debug_flag = 1; -# } - print STDERR "read1 $r1\n" if ($debug_flag == 1); - my $r1sam_flag = $tmp_r1[1]; - my $r2sam_flag = $tmp_r2[1]; - unless ($r1sam_flag) { - print STDERR "error $r1 $r2\n"; - } - next if ($r1sam_flag == 77 || $r1sam_flag == 141); - - my $frag_strand; -### This section is for only properly paired reads - if ($r1sam_flag == 99 || $r1sam_flag == 355) { - $frag_strand = "-"; - } elsif ($r1sam_flag == 83 || $r1sam_flag == 339) { - $frag_strand = "+"; - } elsif ($r1sam_flag == 147 || $r1sam_flag == 403) { - $frag_strand = "-"; - } elsif ($r1sam_flag == 163 || $r1sam_flag == 419) { - $frag_strand = "+"; - } else { - next; - print STDERR "R1 strand error $r1sam_flag\n"; - } -### - # 77 = R1, unmapped - # 141 = R2, unmapped - - # 99 = R1, mapped, fwd strand --- frag on rev strand -> 355 = not primary - # 147 = R2, mapped, rev strand -- frag on rev strand -> 403 = not primary - - # 101 = R1 unmapped, R2 mapped rev strand -- frag on rev strand - # 73 = R1, mapped, fwd strand --- frag on rev strand - # 153 = R2 mapped (R1 unmapped), rev strand -- frag on rev strand - # 133 = R2 unmapped, R1 mapped fwd strand -- frag on rev strand - - # 83 = R1, mapped, rev strand --- frag on fwd strand -> 339 = not primary - # 163 = R2, mapped, fwd strand -- frag on fwd strand -> 419 = not primary - - # 69 = R1 unmapped, R2 mapped fwd strand -- frag on fwd strand - # 89 = R1 mapped rev strand, R2 unmapped -- frag on fwd strand - # 137 = R2 mapped (R1 unmapped), fwd strand -- frag on fwd strand - # 165 = R2 unmapped, R1 rev strand -- frag on fwd strand - - my $flags_r1 = join("\t",@tmp_r1[11..$#tmp_r1]); - my $flags_r2 = join("\t",@tmp_r2[11..$#tmp_r2]); - - my ($mismatch_score_r1,$mismatch_score_r2); - if ($flags_r1 =~ /AS\:i\:(\S+?)\s/) { - $mismatch_score_r1 = $1; - } - if ($flags_r2 =~ /AS\:i\:(\S+?)\s/) { - $mismatch_score_r2 = $1; - } - my $paired_mismatch_score = $mismatch_score_r1 + $mismatch_score_r2; - - - my $mapped_enst = $tmp_r1[2]; - my $mapped_enst_full = $tmp_r1[2]; - if ($mapped_enst =~ /^(.+)\_spliced/) { - $mapped_enst = $1; - } - if ($mapped_enst =~ /^(.+)\_withgenomeflank/) { - $mapped_enst = $1; - } - - unless (exists $convert_enst2type{$mapped_enst}) { - print STDERR "enst2type is missing for $mapped_enst $r1\n"; - } - my ($ensttype,$enstpriority) = split(/\:/,$convert_enst2type{$mapped_enst}); - - if ($frag_strand eq "-") { - $ensttype = "antisense_".$ensttype; - $mapped_enst_full = "antisense_".$mapped_enst_full; - } - - print "mapped $mapped_enst ensttype $ensttype priority $enstpriority\n" if ($debug_flag == 1); - - if ($r1name eq $prev_r1name) { - # current read is same as previous - do nothing - } else { - #current read is different than previous - print and clear hash every 100,000 entries - - $read_counter++; - - if ($read_counter > $print_batch) { - &print_output(); - %read_hash = (); - $read_counter = 0; - } - $prev_r1name = $r1name; - } - - unless (exists $read_hash{$r1name}) { - # if read has never been seen before, keep first mapping - $read_hash{$r1name}{R1}{$ensttype} = $r1; - $read_hash{$r1name}{R2}{$ensttype} = $r2; - $read_hash{$r1name}{flags}{$ensttype} = $enstpriority; - $read_hash{$r1name}{quality} = $paired_mismatch_score; -# push @{$read_hash{$r1name}{mult_ensts}{$ensttype}},$mapped_enst; - push @{$read_hash{$r1name}{mult_ensts}{$ensttype}},$mapped_enst_full; - $read_hash{$r1name}{enst}{$mapped_enst} = $mapped_enst_full; - $read_hash{$r1name}{master_enst}{$ensttype} = $mapped_enst_full; - print "read never seen before quality $paired_mismatch_score $ensttype\n" if ($debug_flag == 1); - - } else { - # is score better than previous mapping? - if ($paired_mismatch_score < $read_hash{$r1name}{quality}) { - # new one is worse, skip - print "new match has worse score than previous - skip $paired_mismatch_score\n" if ($debug_flag == 1); - } elsif ($paired_mismatch_score > $read_hash{$r1name}{quality}) { - # new one is better match than previous - old should all be discarded - delete($read_hash{$r1name}); - $read_hash{$r1name}{R1}{$ensttype} = $r1; - $read_hash{$r1name}{R2}{$ensttype} = $r2; - $read_hash{$r1name}{flags}{$ensttype} = $enstpriority; - $read_hash{$r1name}{quality} = $paired_mismatch_score; -# push @{$read_hash{$r1name}{mult_ensts}{$ensttype}},$mapped_enst; - push @{$read_hash{$r1name}{mult_ensts}{$ensttype}},$mapped_enst_full; - $read_hash{$r1name}{enst}{$mapped_enst} = $mapped_enst_full; - $read_hash{$r1name}{master_enst}{$ensttype} = $mapped_enst_full; - - print "new match has better quality score - discard all previous $paired_mismatch_score\n" if ($debug_flag == 1); - - } elsif ($paired_mismatch_score == $read_hash{$r1name}{quality}) { - # equal quality, both are good - now do family analysis -# print STDERR "two mapping same quality $r1name priority $enstpriority $read_hash{$r1name}{flags}{$ensttype}\n"; - print "has equal quality, now doing family mapping... \n" if ($debug_flag == 1); - - if (exists $read_hash{$r1name}{flags}{$ensttype}) { - # if mapping within family exists before... - print "mapping exists within family before... \n" if ($debug_flag == 1); - # first - did it already map to this transcript before? - if (exists $read_hash{$r1name}{enst}{$mapped_enst}) { - # is it spliced vs unspliced or tRNA flank vs whole genome? if yes ok - if ($read_hash{$r1name}{enst}{$mapped_enst}."_withgenomeflank" eq $mapped_enst_full || $read_hash{$r1name}{enst}{$mapped_enst}."_spliced" eq $mapped_enst_full) { - #original entry is to shorter transcript - keep that one, skip new one entirely - print "original was to original (non-genome flank version for tRNA or non-spliced for others) - keep old, skip this one\n" if ($debug_flag == 1); - } elsif ($read_hash{$r1name}{enst}{$mapped_enst} eq $mapped_enst_full."_withgenomeflank" || $read_hash{$r1name}{enst}{$mapped_enst} eq $mapped_enst_full."_spliced") { - # original entry is to longer transcript - replace with shorter - $read_hash{$r1name}{R1}{$ensttype} = $r1; - $read_hash{$r1name}{R2}{$ensttype} = $r2; - $read_hash{$r1name}{flags}{$ensttype} = $enstpriority; - -# Can't do this as normal - need to replace old -# unshift @{$read_hash{$r1name}{mult_ensts}{$ensttype}},$mapped_enst_full; - unless (exists $read_hash{$r1name}{"mult_ensts"}{$ensttype}) { - print STDERR "error doesn't exist? $r1name $ensttype \n"; - } - for (my $i=0;$i<@{$read_hash{$r1name}{"mult_ensts"}{$ensttype}};$i++) { - if ($read_hash{$r1name}{mult_ensts}{$ensttype}[$i] eq $read_hash{$r1name}{enst}{$mapped_enst}) { - $read_hash{$r1name}{mult_ensts}{$ensttype}[$i] = $mapped_enst_full; - } - } - $read_hash{$r1name}{enst}{$mapped_enst} = $mapped_enst_full; - $read_hash{$r1name}{master_enst}{$ensttype} = $mapped_enst_full; - - print "new mapping is the one i want - discard old and keep newer annotation \n" if ($debug_flag == 1); - } else { - # maps to two places in the same transcript - this is probably actually ok for counting purposes, but for now I'm going to flag these as bad - # 7/20/16 - going to comment this out for now - basically just skip the second entry but don't flag as bad -# $read_hash{$r1name}{flags}{"double_maps"} = 1; - - print "double maps - $ensttype prev length of mult_ensts array is ".scalar(@{$read_hash{$r1name}{mult_ensts}{$ensttype}})."\n" if ($debug_flag == 1); - - for (my $i=0;$i<@{$read_hash{$r1name}{"mult_ensts"}{$ensttype}};$i++) { - if ($read_hash{$r1name}{mult_ensts}{$ensttype}[$i] eq $read_hash{$r1name}{enst}{$mapped_enst}) { - $read_hash{$r1name}{mult_ensts}{$ensttype}[$i] = $mapped_enst_full."_DOUBLEMAP"; - } - } - - $read_hash{$r1name}{master_enst}{$ensttype} = $mapped_enst_full."_DOUBLEMAP"; - - print "double maps - $ensttype length of mult_ensts array is ".scalar(@{$read_hash{$r1name}{mult_ensts}{$ensttype}})."\n" if ($debug_flag == 1); - } - - - } elsif ($enstpriority < $read_hash{$r1name}{flags}{$ensttype}) { - # priority of new mapping is better than old - replace old - print "priority of new mapping is better than old - replace old mapping\n" if ($debug_flag == 1); - $read_hash{$r1name}{R1}{$ensttype} = $r1; - $read_hash{$r1name}{R2}{$ensttype} = $r2; - $read_hash{$r1name}{flags}{$ensttype} = $enstpriority; -# unshift @{$read_hash{$r1name}{mult_ensts}{$ensttype}},$mapped_enst; - unshift @{$read_hash{$r1name}{mult_ensts}{$ensttype}},$mapped_enst_full; -# print STDERR "adding mapped $r1name $mapped_enst\n"; - - $read_hash{$r1name}{enst}{$mapped_enst} = $mapped_enst_full; - $read_hash{$r1name}{master_enst}{$ensttype} = $mapped_enst_full; - } else { - # Mapping is equal quality, but priority of new mapping is worse than old - keep new enst_full but otherwise discard -# push @{$read_hash{$r1name}{mult_ensts}{$ensttype}},$mapped_enst; - push @{$read_hash{$r1name}{mult_ensts}{$ensttype}},$mapped_enst_full; - } - } elsif (exists $rRNA_extra_hash{$ensttype} && exists $read_hash{$r1name}{R1}{$rRNA_extra_hash{$ensttype}}) { - - # If old mapping was to rRNA_extra but new mapping is RNA18S, RNA28S, or RNA5-8S, keep new mapping and discard old - my $old_rRNA_label = $rRNA_extra_hash{$ensttype}; - delete($read_hash{$r1name}{R1}{$old_rRNA_label}); - delete($read_hash{$r1name}{R2}{$old_rRNA_label}); - delete($read_hash{$r1name}{flags}{$old_rRNA_label}); - delete($read_hash{$r1name}{master_enst}{$old_rRNA_label}); - delete($read_hash{$r1name}{mult_ensts}{$old_rRNA_label}); - delete($read_hash{$r1name}{enst}{"NR_046235.1"}); - - $read_hash{$r1name}{R1}{$ensttype} = $r1; - $read_hash{$r1name}{R2}{$ensttype} = $r2; - $read_hash{$r1name}{flags}{$ensttype} = $enstpriority; - $read_hash{$r1name}{quality} = $paired_mismatch_score; - $read_hash{$r1name}{enst}{$mapped_enst} = $mapped_enst_full; - $read_hash{$r1name}{master_enst}{$ensttype} = $mapped_enst_full; - push @{$read_hash{$r1name}{mult_ensts}{$ensttype}},$mapped_enst_full; - - print "new mapping is the one i want - discard old and keep newer annotation \n" if ($debug_flag == 1); - } elsif (exists $rRNA_extra_hash_rev{$ensttype}) { - # new mapping is to rRNA_extra - my $rRNA_flag = 0; - for my $rRNA_elements (keys %{$rRNA_extra_hash_rev{$ensttype}}) { - if (exists $read_hash{$r1name}{R1}{$rRNA_elements}) { - # new mapping is to rRNA_extra, but old is to RNA28S, RNA18S, or RNA5-8S - discard new mapping - $rRNA_flag = 1; - } - } - - if ($rRNA_flag == 1) {} else { - # new mapping is to rRNA_extra, but old is to something non-rRNA - do same as multi-family below - $read_hash{$r1name}{R1}{$ensttype} = $r1; - $read_hash{$r1name}{R2}{$ensttype} = $r2; - $read_hash{$r1name}{flags}{$ensttype} = $enstpriority; - $read_hash{$r1name}{quality} = $paired_mismatch_score; - $read_hash{$r1name}{enst}{$mapped_enst} = $mapped_enst_full; - push @{$read_hash{$r1name}{mult_ensts}{$ensttype}},$mapped_enst_full; - $read_hash{$r1name}{master_enst}{$ensttype} = $mapped_enst_full; - - } -# $rRNA_extra_hash_rev{"rRNA_extra"}{"RNA28S"} = 1; - - } else { - print "maps to two families - $ensttype\n" if ($debug_flag == 1); -# hits this if - # 1) read_hash{r1name} already exists (read already seen) - # 2) $paired_mismatch_score == $read_hash{$r1name}{quality} - mapping score is same as previous - # 3) $read_hash{$r1name}{flags}{$ensttype} doesn't exist - read exists but this is a new family - - # read maps to multiple gene lists - discard as not uniquely mapping within family - # 7/19/16 - changed this to KEEP multi-family mapping - - - #my $type1 = "RN7SL"; - #my $type2 = "RN7SL_antisense"; -# if (exists $read_hash{$r1name}{flags}{$type1} && $ensttype eq $type2) { -# print "$type1\t".$read_hash{$r1name}{R1}{$type1}."\n".$read_hash{$r1name}{R2}{$type1}."\n"."$type2\t".$r1."\n".$r2."\n\n"; -# } -# if (exists $read_hash{$r1name}{flags}{$type2} && $ensttype eq $type1) { -# print "$type2\t".$read_hash{$r1name}{R1}{$type2}."\n".$read_hash{$r1name}{R2}{$type2}."\n"."$type1\t".$r1."\n".$r2."\n\n"; -# } - - $read_hash{$r1name}{R1}{$ensttype} = $r1; - $read_hash{$r1name}{R2}{$ensttype} = $r2; - $read_hash{$r1name}{flags}{$ensttype} = $enstpriority; - $read_hash{$r1name}{quality} = $paired_mismatch_score; - $read_hash{$r1name}{enst}{$mapped_enst} = $mapped_enst_full; - push @{$read_hash{$r1name}{mult_ensts}{$ensttype}},$mapped_enst_full; - - #adds new master enst to save for below - $read_hash{$r1name}{master_enst}{$ensttype} = $mapped_enst_full; -# my @tmp_to_sort = split(/\|/,$read_hash{$r1name}{master_enst}); -# push @tmp_to_sort,$mapped_enst_full; -# my @sorted = sort {$a cmp $b} @tmp_to_sort; -# $read_hash{$r1name}{master_enst} = join("|",@sorted); - - } - } else { - # I don't think this should ever be hit - print STDERR "this shouldn't be hit - $paired_mismatch_score $read_hash{$r1name}{quality} $r1name $r1 $r2\n"; - } - - } - # end of line processing - - } -} - -&print_output(); -close(SAMOUT); -for my $multi_key (keys %multimapping_hash) { - print MULTIMAP "$multi_key\t".$multimapping_hash{$multi_key}."\n"; -} -close(MULTIMAP); - -#my $done_file = $output.".done"; -open(DONE,">$done_file"); -print DONE "jobs done\n"; -close(DONE); - - -sub print_output { - my %count; - my %count_enst; - for my $read (keys %read_hash) { - my @ensttype_array = sort {$a cmp $b} keys %{$read_hash{$read}{flags}}; - my $ensttype = $ensttype_array[0]; - my @masterenst_array; - for my $type (@ensttype_array) { - push @masterenst_array,$read_hash{$read}{master_enst}{$type}; - } - my $ensttype_join = join("|",@ensttype_array); - my $masterenst_join = join("|",@masterenst_array); - - $count{$ensttype_join}++; - -# print STDERR "read $read, ensttype $ensttype\n"; - if (scalar(keys %{$read_hash{$read}{flags}}) == 1) { - print STDERR "this shouldn't happen this should be 1 ".scalar(keys %{$read_hash{$read}{R1}})."\n" unless (scalar(keys %{$read_hash{$read}{R1}}) == 1); - - my @r1_split = split(/\t/,$read_hash{$read}{R1}{$ensttype}); - my @r2_split = split(/\t/,$read_hash{$read}{R2}{$ensttype}); - $r1_split[2] = $ensttype_join."||".$masterenst_join; - $r2_split[2] = $ensttype_join."||".$masterenst_join; - my $r1_line = join("\t",@r1_split); - my $r2_line = join("\t",@r2_split); - - print SAMOUT "".$r1_line."\tZZ:Z:".join("|",@{$read_hash{$read}{mult_ensts}{$ensttype}})."\n".$r2_line."\tZZ:Z:".join("|",@{$read_hash{$read}{mult_ensts}{$ensttype}})."\n"; -# print SAMOUT "".$read_hash{$read}{R1}{$ensttype}."\tZZ:Z:".join("|",@{$read_hash{$read}{mult_ensts}{$ensttype}})."\n".$read_hash{$read}{R2}{$ensttype}."\tZZ:Z:".join("|",@{$read_hash{$read}{mult_ensts}{$ensttype}})."\n"; - - my @blah = split(/\t/,$read_hash{$read}{R1}{$ensttype}); - $count_enst{$ensttype."|".$blah[2]}++; -# $count_enst{$ensttype."|".join("|",@{$read_hash{$read}{mult_ensts}{$ensttype}})}++; - - } else { -# print STDERR "this shouldn't happen this should be 1 ".scalar(keys %{$read_hash{$read}{R1}})."\n" unless (scalar(keys %{$read_hash{$read}{R1}}) == 1); - my @all_mult_ensts; - for my $key (@ensttype_array) { - push @all_mult_ensts,join("|",@{$read_hash{$read}{mult_ensts}{$key}}); - } - my $final_mult_ensts = join("|",@all_mult_ensts); - - unless (exists $read_hash{$read}{R1}{$ensttype} && $read_hash{$read}{R1}{$ensttype}) { - print "weird error - $read $ensttype readhash doesn't exist ? ".$read_hash{$read}{flags}{$ensttype}."\n"; - } - - my @r1_split = split(/\t/,$read_hash{$read}{R1}{$ensttype}); - my @r2_split = split(/\t/,$read_hash{$read}{R2}{$ensttype}); - $r1_split[2] = $ensttype_join."||".$masterenst_join; - $r2_split[2] = $ensttype_join."||".$masterenst_join; - - my $r1_line = join("\t",@r1_split); - my $r2_line = join("\t",@r2_split); -# print SAMOUT "".$r1_line."\tZZ:Z:".join("|",@{$read_hash{$read}{mult_ensts}{$ensttype}})."\n".$r2_line."\tZZ:Z:".join("|",@{$read_hash{$read}{mult_ensts}{$ensttype}})."\n"; - print SAMOUT "".$r1_line."\tZZ:Z:".$final_mult_ensts."\n".$r2_line."\tZZ:Z:".$final_mult_ensts."\n"; - my $multimapping_type = join("|",keys %{$read_hash{$read}{flags}}); - $multimapping_hash{$multimapping_type}++; - } - } - - -} - - -#my @sorted = sort {$count{$b} <=> $count{$a}} keys %count; -#for my $k (@sorted) { -# print "$k $count{$k}\n"; -#} - -#my @sorted = sort {$count_enst{$b} <=> $count_enst{$a}} keys %count_enst; -#for my $s (@sorted) { -# my ($type,$geneid) = split(/\|/,$s); -# print "$type\t$geneid\t$count_enst{$s}\t$enst2gene{$geneid}\n"; -# print "$type\t$geneid\t$count_enst{$s}\t$s\n"; -#} - - -#print "unique $unique_count\n duplicate $duplicate_count\n"; - - - - - -sub read_in_filelists { - my $fi = shift; - my $priority_n = 0; -# print STDERR "opening $fi\n"; - open(F,$fi); - for my $line () { - chomp($line); - my ($allenst,$allensg,$gid,$type_label,$typefile) = split(/\t/,$line); - unless ($type_label) { - print "err $line\n"; - } -# my ($allensg,$gid,$allenst) = split(/\t/,$line); - $type_label =~ s/\_$//; - unless ($allenst) { - print STDERR "error missing enst $line $fi\n"; - } - my @ensts = split(/\|/,$allenst); - for my $enst (@ensts) { - $enst2gene{$enst} = $gid; - $convert_enst2type{$enst} = $type_label.":".$priority_n; - $priority_n++; - } - } - close(F); -} diff --git a/bin/perl/parse_bowtie2_output_realtime_includemultifamily_SE.pl b/bin/perl/parse_bowtie2_output_realtime_includemultifamily_SE.pl deleted file mode 100755 index be9c36e..0000000 --- a/bin/perl/parse_bowtie2_output_realtime_includemultifamily_SE.pl +++ /dev/null @@ -1,481 +0,0 @@ -#!/usr/bin/env perl - -use warnings; -use strict; - -my %enst2gene; -my %convert_enst2type; -my %multimapping_hash; - -my %rRNA_extra_hash; -my %rRNA_extra_hash_rev; -$rRNA_extra_hash{"RNA28S"} = "RNA45S"; -$rRNA_extra_hash{"RNA18S"} = "RNA45S"; -$rRNA_extra_hash{"RNA5-8S"} = "RNA45S"; -$rRNA_extra_hash{"antisense_RNA28S"} = "antisense_RNA45S"; -$rRNA_extra_hash{"antisense_RNA18S"} = "antisense_RNA45S"; -$rRNA_extra_hash{"antisense_RNA5-8S"} = "antisense_RNA45S"; - -$rRNA_extra_hash_rev{"RNA45S"}{"RNA28S"} = 1; -$rRNA_extra_hash_rev{"RNA45S"}{"RNA18S"} = 1; -$rRNA_extra_hash_rev{"RNA45S"}{"RNA5-8S"} = 1; -$rRNA_extra_hash_rev{"antisense_RNA45S"}{"antisense_RNA28S"} = 1; -$rRNA_extra_hash_rev{"antisense_RNA45S"}{"antisense_RNA18S"} = 1; -$rRNA_extra_hash_rev{"antisense_RNA45S"}{"antisense_RNA5-8S"} = 1; - - - -# my $species = $ARGV[3]; -#my $working_directory = "/home/elvannostrand/data/clip/CLIPseq_analysis/ENCODE_v9_20151209/PolIII_mapping/"; -# my $working_directory = $ARGV[2]; - -my $filelist_file = $ARGV[3]; -# my $filelist_file = "/home/elvannostrand/data/clip/CLIPseq_analysis/RNA_type_analysis/MASTER_filelist.wrepbaseandtRNA.enst2id.fixed.UpdatedSimpleRepeat.20190424"; -# if ($species eq "hg38") { -# $filelist_file = "/home/elvannostrand/data/clip/CLIPseq_analysis/RNA_type_analysis/hg38/MASTER_FILELIST.20201203.wrepbaseandtRNA.enst2id.fixed.UpdatedSimpleRepeat.wmiRs.list"; -# } - -&read_in_filelists($filelist_file); - -my $print_batch = 10000; -my $read_counter = 0; - -#&read_in_filelists("/home/elvannostrand/data/clip/CLIPseq_analysis/RNA_type_analysis/genelists.RNA5S","RNA5S"); -#my $fastq_file1 = "/home/gpratt/projects/encode/analysis/encode_v9/484_CLIP_S6_L001_R1_001.B06_484_01_POLR2G.adapterTrim.round2.fastq.gz"; -#my $fastq_file2 = "/home/gpratt/projects/encode/analysis/encode_v9/484_CLIP_S6_L001_R2_001.B06_484_01_POLR2G.adapterTrim.round2.fastq.gz"; - -#my $fastq_file1 = "test_r1.fastq"; -#my $fastq_file2 = "test_r2.fastq"; -#my $fastq_file1 = "484_CLIP_S6_L001_R1_001.B06_484_01_POLR2G.adapterTrim.round2.fastq"; -#my $fastq_file2 = "484_CLIP_S6_L001_R2_001.B06_484_01_POLR2G.adapterTrim.round2.fastq"; -#my $fastq_file1 = "484_INPUT_S7_L001_R1_001.unassigned.adapterTrim.round2.fastq"; -#my $fastq_file2 = "484_INPUT_S7_L001_R2_001.unassigned.adapterTrim.round2.fastq"; - -my $fastq_file1 = $ARGV[0]; - -my $bowtie_db = $ARGV[1]; -#my $bowtie_db = "/home/elvannostrand/data/clip/CLIPseq_analysis/RNA_type_analysis/filelist_POLIII.combined.Nflank"; -my @bowtie_db_split = split(/\//,$bowtie_db); -my $bowtie_db_short = $bowtie_db_split[$#bowtie_db_split]; - -my @fastq_fi1_split = split(/\//,$fastq_file1); -my $fastq_fi1_short = $fastq_fi1_split[$#fastq_fi1_split]; - -my $output = $ARGV[2]; -#my $output = $working_directory.$fastq_fi1_short.".mapped_vs_".$bowtie_db_short.".sam"; -my $bowtie_out = $output.".bowtieout"; -open(SAMOUT,">$output"); -my $multimapping_out = $output.".multimapping_deleted"; -open(MULTIMAP,">$multimapping_out"); -my $done_file = $output.".done"; - -# -my $command = "stdbuf -oL bowtie2 -q --sensitive -a -p 1 --no-mixed --reorder -x $bowtie_db -U $fastq_file1 2> $bowtie_out"; -#my $command = "/projects/ps-yeolab/software/bowtie-1.1.1/./bowtie -q -1 $fastq_file1 -2 $fastq_file2 $bowtie_db -a -v 2 --best --strata -S 2> $bowtie_out"; -print STDERR "command $command\n"; -#my $pid = open(BOWTIE, "-|", "unbuffer /projects/ps-yeolab/software/bowtie-1.1.1/./bowtie -q -1 $fastq_file1 -2 $fastq_file2 $genelistfi -a -v 2 --best --strata -S 2> $bowtie_out"); -my $pid = open(BOWTIE, "-|", "stdbuf -oL bowtie2 -q --sensitive -a -p 3 --no-mixed --reorder -x $bowtie_db -U $fastq_file1 2> $bowtie_out"); -#my $pid = open(BOWTIE, "-|", "stdbuf -oL /projects/ps-yeolab/software/bowtie-1.1.1/./bowtie -q -1 $fastq_file1 -2 $fastq_file2 $bowtie_db -a -v 2 --best --strata -S 2> $bowtie_out"); -#print "PID $pid\n"; -my %fragment_hash; -my $duplicate_count=0; -my $unique_count=0; -my $all_count=0; -my $prev_r1name = ""; - -my %read_hash; -if ($pid) { - - while () { - my $r1 = $_; - chomp($r1); - if ($r1 =~ /^\@/) { - print SAMOUT "$r1\n"; - next; - } - - my @tmp_r1 = split(/\t/,$r1); - -# my ($r1name,$r1bc) = split(/\s+/,$tmp_r1[0]); -# my ($r1name,$r1bc) = split(/\_/,$tmp_r1[0]); - my @read_name = split(/\_/,$tmp_r1[0]); - my $r1bc = pop(@read_name); - my $r1name = join("_",@read_name); - -# $all_count++; -# if ($all_count % 1000 == 0) { -# print STDERR "$all_count\n"; -# } - - my $debug_flag = 0; -# if ($r1name eq "K00180:734:H2FYCBBXY:5:1101:12053:6378_TCCCCCCGAC") { -# print "READ1\t$r1\n"; -# $debug_flag = 1; -# } - print STDERR "read1 $r1\n" if ($debug_flag == 1); - my $r1sam_flag = $tmp_r1[1]; - -# unless ($r1sam_flag) { -# print STDERR "error $r1\n"; -# } - next if ($r1sam_flag == 4); - - my $frag_strand; -### This section is for only properly paired reads - if ($r1sam_flag == 16 || $r1sam_flag == 272) { - $frag_strand = "-"; - } elsif ($r1sam_flag eq "0" || $r1sam_flag == 256) { - $frag_strand = "+"; - } else { - next; - print STDERR "R1 strand error $r1sam_flag\n"; - } -### - my $flags_r1 = join("\t",@tmp_r1[11..$#tmp_r1]); - - my ($mismatch_score_r1); - if ($flags_r1 =~ /AS\:i\:(\S+?)\s/) { - $mismatch_score_r1 = $1; - } - my $paired_mismatch_score = $mismatch_score_r1; - - - my $mapped_enst = $tmp_r1[2]; - my $mapped_enst_full = $tmp_r1[2]; - if ($mapped_enst =~ /^(.+)\_spliced/) { - $mapped_enst = $1; - } - if ($mapped_enst =~ /^(.+)\_withgenomeflank/) { - $mapped_enst = $1; - } - - unless (exists $convert_enst2type{$mapped_enst}) { - print STDERR "enst2type is missing for $mapped_enst $r1\n"; - } - my ($ensttype,$enstpriority) = split(/\:/,$convert_enst2type{$mapped_enst}); - - if ($frag_strand eq "-") { - $ensttype = "antisense_".$ensttype; - $mapped_enst_full = "antisense_".$mapped_enst_full; - } - - print "mapped $mapped_enst ensttype $ensttype priority $enstpriority\n" if ($debug_flag == 1); - - if ($r1name eq $prev_r1name) { - # current read is same as previous - do nothing - } else { - #current read is different than previous - print and clear hash every 100,000 entries - - $read_counter++; - - if ($read_counter > $print_batch) { - &print_output(); - %read_hash = (); - $read_counter = 0; - } - $prev_r1name = $r1name; - } - - if ($debug_flag == 1) { - print "r1name $r1name \t paired_mismatch_score $paired_mismatch_score\n"; - if (exists $read_hash{$r1name}) { - print "read has been seen before - $read_hash{$r1name}{quality}\n"; - } else { - print "read never seen before\n"; - } - } - - unless (exists $read_hash{$r1name}) { - # if read has never been seen before, keep first mapping - $read_hash{$r1name}{R1}{$ensttype} = $r1; - $read_hash{$r1name}{flags}{$ensttype} = $enstpriority; - $read_hash{$r1name}{quality} = $paired_mismatch_score; -# push @{$read_hash{$r1name}{mult_ensts}{$ensttype}},$mapped_enst; - push @{$read_hash{$r1name}{mult_ensts}{$ensttype}},$mapped_enst_full; - $read_hash{$r1name}{enst}{$mapped_enst} = $mapped_enst_full; - $read_hash{$r1name}{master_enst}{$ensttype} = $mapped_enst_full; - print "read never seen before quality $paired_mismatch_score $ensttype\n" if ($debug_flag == 1); - - } else { - # is score better than previous mapping? - if ($paired_mismatch_score < $read_hash{$r1name}{quality}) { - # new one is worse, skip - print "new match has worse score than previous - skip $paired_mismatch_score\n" if ($debug_flag == 1); - } elsif ($paired_mismatch_score > $read_hash{$r1name}{quality}) { - # new one is better match than previous - old should all be discarded - delete($read_hash{$r1name}); - $read_hash{$r1name}{R1}{$ensttype} = $r1; - $read_hash{$r1name}{flags}{$ensttype} = $enstpriority; - $read_hash{$r1name}{quality} = $paired_mismatch_score; -# push @{$read_hash{$r1name}{mult_ensts}{$ensttype}},$mapped_enst; - push @{$read_hash{$r1name}{mult_ensts}{$ensttype}},$mapped_enst_full; - $read_hash{$r1name}{enst}{$mapped_enst} = $mapped_enst_full; - $read_hash{$r1name}{master_enst}{$ensttype} = $mapped_enst_full; - - print "new match has better quality score - discard all previous $paired_mismatch_score\n" if ($debug_flag == 1); - - } elsif ($paired_mismatch_score == $read_hash{$r1name}{quality}) { - # equal quality, both are good - now do family analysis -# print STDERR "two mapping same quality $r1name priority $enstpriority $read_hash{$r1name}{flags}{$ensttype}\n"; - print "has equal quality, now doing family mapping... \n" if ($debug_flag == 1); - - if (exists $read_hash{$r1name}{flags}{$ensttype}) { - # if mapping within family exists before... - print "mapping exists within family before... \n" if ($debug_flag == 1); - # first - did it already map to this transcript before? - if (exists $read_hash{$r1name}{enst}{$mapped_enst}) { - # is it spliced vs unspliced or tRNA flank vs whole genome? if yes ok - if ($read_hash{$r1name}{enst}{$mapped_enst}."_withgenomeflank" eq $mapped_enst_full || $read_hash{$r1name}{enst}{$mapped_enst}."_spliced" eq $mapped_enst_full) { - #original entry is to shorter transcript - keep that one, skip new one entirely - print "original was to original (non-genome flank version for tRNA or non-spliced for others) - keep old, skip this one\n" if ($debug_flag == 1); - } elsif ($read_hash{$r1name}{enst}{$mapped_enst} eq $mapped_enst_full."_withgenomeflank" || $read_hash{$r1name}{enst}{$mapped_enst} eq $mapped_enst_full."_spliced") { - # original entry is to longer transcript - replace with shorter - $read_hash{$r1name}{R1}{$ensttype} = $r1; - $read_hash{$r1name}{flags}{$ensttype} = $enstpriority; - -# Can't do this as normal - need to replace old -# unshift @{$read_hash{$r1name}{mult_ensts}{$ensttype}},$mapped_enst_full; - unless (exists $read_hash{$r1name}{"mult_ensts"}{$ensttype}) { - print STDERR "error doesn't exist? $r1name $ensttype \n"; - } - for (my $i=0;$i<@{$read_hash{$r1name}{"mult_ensts"}{$ensttype}};$i++) { - if ($read_hash{$r1name}{mult_ensts}{$ensttype}[$i] eq $read_hash{$r1name}{enst}{$mapped_enst}) { - $read_hash{$r1name}{mult_ensts}{$ensttype}[$i] = $mapped_enst_full; - } - } - $read_hash{$r1name}{enst}{$mapped_enst} = $mapped_enst_full; - $read_hash{$r1name}{master_enst}{$ensttype} = $mapped_enst_full; - - print "new mapping is the one i want - discard old and keep newer annotation \n" if ($debug_flag == 1); - } else { - # maps to two places in the same transcript - this is probably actually ok for counting purposes, but for now I'm going to flag these as bad - # 7/20/16 - going to comment this out for now - basically just skip the second entry but don't flag as bad -# $read_hash{$r1name}{flags}{"double_maps"} = 1; - - print "double maps - $ensttype prev length of mult_ensts array is ".scalar(@{$read_hash{$r1name}{mult_ensts}{$ensttype}})."\n" if ($debug_flag == 1); - - for (my $i=0;$i<@{$read_hash{$r1name}{"mult_ensts"}{$ensttype}};$i++) { - if ($read_hash{$r1name}{mult_ensts}{$ensttype}[$i] eq $read_hash{$r1name}{enst}{$mapped_enst}) { - $read_hash{$r1name}{mult_ensts}{$ensttype}[$i] = $mapped_enst_full."_DOUBLEMAP"; - } - } - - $read_hash{$r1name}{master_enst}{$ensttype} = $mapped_enst_full."_DOUBLEMAP"; - - print "double maps - $ensttype length of mult_ensts array is ".scalar(@{$read_hash{$r1name}{mult_ensts}{$ensttype}})."\n" if ($debug_flag == 1); - } - - - } elsif ($enstpriority < $read_hash{$r1name}{flags}{$ensttype}) { - # priority of new mapping is better than old - replace old - print "priority of new mapping is better than old - replace old mapping\n" if ($debug_flag == 1); - $read_hash{$r1name}{R1}{$ensttype} = $r1; - $read_hash{$r1name}{flags}{$ensttype} = $enstpriority; -# unshift @{$read_hash{$r1name}{mult_ensts}{$ensttype}},$mapped_enst; - unshift @{$read_hash{$r1name}{mult_ensts}{$ensttype}},$mapped_enst_full; -# print STDERR "adding mapped $r1name $mapped_enst\n"; - - $read_hash{$r1name}{enst}{$mapped_enst} = $mapped_enst_full; - $read_hash{$r1name}{master_enst}{$ensttype} = $mapped_enst_full; - } else { - # Mapping is equal quality, but priority of new mapping is worse than old - keep new enst_full but otherwise discard -# push @{$read_hash{$r1name}{mult_ensts}{$ensttype}},$mapped_enst; - push @{$read_hash{$r1name}{mult_ensts}{$ensttype}},$mapped_enst_full; - } - } elsif (exists $rRNA_extra_hash{$ensttype} && exists $read_hash{$r1name}{R1}{$rRNA_extra_hash{$ensttype}}) { - - # If old mapping was to rRNA_extra but new mapping is RNA18S, RNA28S, or RNA5-8S, keep new mapping and discard old - my $old_rRNA_label = $rRNA_extra_hash{$ensttype}; - delete($read_hash{$r1name}{R1}{$old_rRNA_label}); - delete($read_hash{$r1name}{flags}{$old_rRNA_label}); - delete($read_hash{$r1name}{master_enst}{$old_rRNA_label}); - delete($read_hash{$r1name}{mult_ensts}{$old_rRNA_label}); - delete($read_hash{$r1name}{enst}{"NR_046235.1"}); - - $read_hash{$r1name}{R1}{$ensttype} = $r1; - $read_hash{$r1name}{flags}{$ensttype} = $enstpriority; - $read_hash{$r1name}{quality} = $paired_mismatch_score; - $read_hash{$r1name}{enst}{$mapped_enst} = $mapped_enst_full; - $read_hash{$r1name}{master_enst}{$ensttype} = $mapped_enst_full; - push @{$read_hash{$r1name}{mult_ensts}{$ensttype}},$mapped_enst_full; - - print "new mapping is the one i want - discard old and keep newer annotation \n" if ($debug_flag == 1); - } elsif (exists $rRNA_extra_hash_rev{$ensttype}) { - # new mapping is to rRNA_extra - my $rRNA_flag = 0; - for my $rRNA_elements (keys %{$rRNA_extra_hash_rev{$ensttype}}) { - if (exists $read_hash{$r1name}{R1}{$rRNA_elements}) { - # new mapping is to rRNA_extra, but old is to RNA28S, RNA18S, or RNA5-8S - discard new mapping - $rRNA_flag = 1; - } - } - - if ($rRNA_flag == 1) {} else { - # new mapping is to rRNA_extra, but old is to something non-rRNA - do same as multi-family below - print STDERR "new mapping is to rRNA_extra, but old is to something non-rRNA - do same as multi-family below\nrRNA_flag $rRNA_flag ensttype $ensttype r1 $r1\n" if ($debug_flag == 1); - $read_hash{$r1name}{R1}{$ensttype} = $r1; - $read_hash{$r1name}{flags}{$ensttype} = $enstpriority; - $read_hash{$r1name}{quality} = $paired_mismatch_score; - $read_hash{$r1name}{enst}{$mapped_enst} = $mapped_enst_full; - push @{$read_hash{$r1name}{mult_ensts}{$ensttype}},$mapped_enst_full; - $read_hash{$r1name}{master_enst}{$ensttype} = $mapped_enst_full; - - } - - - } else { - print "maps to two families - $ensttype\n" if ($debug_flag == 1); -# hits this if - # 1) read_hash{r1name} already exists (read already seen) - # 2) $paired_mismatch_score == $read_hash{$r1name}{quality} - mapping score is same as previous - # 3) $read_hash{$r1name}{flags}{$ensttype} doesn't exist - read exists but this is a new family - - # read maps to multiple gene lists - discard as not uniquely mapping within family - # 7/19/16 - changed this to KEEP multi-family mapping - - - #my $type1 = "RN7SL"; - #my $type2 = "RN7SL_antisense"; -# if (exists $read_hash{$r1name}{flags}{$type1} && $ensttype eq $type2) { -# print "$type1\t".$read_hash{$r1name}{R1}{$type1}."\n".$read_hash{$r1name}{R2}{$type1}."\n"."$type2\t".$r1."\n".$r2."\n\n"; -# } -# if (exists $read_hash{$r1name}{flags}{$type2} && $ensttype eq $type1) { -# print "$type2\t".$read_hash{$r1name}{R1}{$type2}."\n".$read_hash{$r1name}{R2}{$type2}."\n"."$type1\t".$r1."\n".$r2."\n\n"; -# } - - $read_hash{$r1name}{R1}{$ensttype} = $r1; - $read_hash{$r1name}{flags}{$ensttype} = $enstpriority; - $read_hash{$r1name}{quality} = $paired_mismatch_score; - $read_hash{$r1name}{enst}{$mapped_enst} = $mapped_enst_full; - push @{$read_hash{$r1name}{mult_ensts}{$ensttype}},$mapped_enst_full; - - #adds new master enst to save for below - $read_hash{$r1name}{master_enst}{$ensttype} = $mapped_enst_full; -# my @tmp_to_sort = split(/\|/,$read_hash{$r1name}{master_enst}); -# push @tmp_to_sort,$mapped_enst_full; -# my @sorted = sort {$a cmp $b} @tmp_to_sort; -# $read_hash{$r1name}{master_enst} = join("|",@sorted); - - } - } else { - # I don't think this should ever be hit -# print STDERR "this shouldn't be hit - $paired_mismatch_score $read_hash{$r1name}{quality} $r1name $r1 $r2\n"; - print STDERR "this shouldn't be hit - $paired_mismatch_score $read_hash{$r1name}{quality} $r1name $r1\n"; - } - - } - # end of line processing - - } -} - -&print_output(); -close(SAMOUT); -for my $multi_key (keys %multimapping_hash) { - print MULTIMAP "$multi_key\t".$multimapping_hash{$multi_key}."\n"; -} -close(MULTIMAP); - -#my $done_file = $output.".done"; -open(DONE,">$done_file"); -print DONE "jobs done\n"; -close(DONE); - - -sub print_output { - my %count; - my %count_enst; - for my $read (keys %read_hash) { - my @ensttype_array = sort {$a cmp $b} keys %{$read_hash{$read}{flags}}; - my $ensttype = $ensttype_array[0]; - my @masterenst_array; - for my $type (@ensttype_array) { - push @masterenst_array,$read_hash{$read}{master_enst}{$type}; - } - my $ensttype_join = join("|",@ensttype_array); - my $masterenst_join = join("|",@masterenst_array); - - $count{$ensttype_join}++; - -# print STDERR "read $read, ensttype $ensttype\n"; - if (scalar(keys %{$read_hash{$read}{flags}}) == 1) { - print STDERR "this shouldn't happen this should be 1 ".scalar(keys %{$read_hash{$read}{R1}})."\n" unless (scalar(keys %{$read_hash{$read}{R1}}) == 1); - - my @r1_split = split(/\t/,$read_hash{$read}{R1}{$ensttype}); - $r1_split[2] = $ensttype_join."||".$masterenst_join; - my $r1_line = join("\t",@r1_split); - -# print SAMOUT "".$r1_line."\tZZ:Z:".join("|",@{$read_hash{$read}{mult_ensts}{$ensttype}})."\n".$r2_line."\tZZ:Z:".join("|",@{$read_hash{$read}{mult_ensts}{$ensttype}})."\n"; - print SAMOUT "".$r1_line."\tZZ:Z:".join("|",@{$read_hash{$read}{mult_ensts}{$ensttype}})."\n"; - - my @blah = split(/\t/,$read_hash{$read}{R1}{$ensttype}); - $count_enst{$ensttype."|".$blah[2]}++; -# $count_enst{$ensttype."|".join("|",@{$read_hash{$read}{mult_ensts}{$ensttype}})}++; - - } else { -# print STDERR "this shouldn't happen this should be 1 ".scalar(keys %{$read_hash{$read}{R1}})."\n" unless (scalar(keys %{$read_hash{$read}{R1}}) == 1); - my @all_mult_ensts; - for my $key (@ensttype_array) { - push @all_mult_ensts,join("|",@{$read_hash{$read}{mult_ensts}{$key}}); - } - my $final_mult_ensts = join("|",@all_mult_ensts); - - unless (exists $read_hash{$read}{R1}{$ensttype} && $read_hash{$read}{R1}{$ensttype}) { - print "weird error - $read $ensttype readhash doesn't exist ? ".$read_hash{$read}{flags}{$ensttype}."\n"; - } - - my @r1_split = split(/\t/,$read_hash{$read}{R1}{$ensttype}); - $r1_split[2] = $ensttype_join."||".$masterenst_join; - my $r1_line = join("\t",@r1_split); - - print SAMOUT "".$r1_line."\tZZ:Z:".$final_mult_ensts."\n"; -# print SAMOUT "".$r1_line."\tZZ:Z:".$final_mult_ensts."\n".$r2_line."\tZZ:Z:".$final_mult_ensts."\n"; - my $multimapping_type = join("|",keys %{$read_hash{$read}{flags}}); - $multimapping_hash{$multimapping_type}++; - } - } - - -} - - -#my @sorted = sort {$count{$b} <=> $count{$a}} keys %count; -#for my $k (@sorted) { -# print "$k $count{$k}\n"; -#} - -#my @sorted = sort {$count_enst{$b} <=> $count_enst{$a}} keys %count_enst; -#for my $s (@sorted) { -# my ($type,$geneid) = split(/\|/,$s); -# print "$type\t$geneid\t$count_enst{$s}\t$enst2gene{$geneid}\n"; -# print "$type\t$geneid\t$count_enst{$s}\t$s\n"; -#} - - -#print "unique $unique_count\n duplicate $duplicate_count\n"; - - - - - -sub read_in_filelists { - my $fi = shift; - my $priority_n = 0; - open(F,$fi) || die "no $fi\n"; - for my $line () { - chomp($line); - my ($allenst,$allensg,$gid,$type_label,$typefile) = split(/\t/,$line); -# my ($allensg,$gid,$allenst) = split(/\t/,$line); - $type_label =~ s/\_$//; - unless ($allenst) { - print STDERR "error missing enst $line $fi\n"; - } - my @ensts = split(/\|/,$allenst); - for my $enst (@ensts) { - $enst2gene{$enst} = $gid; - $convert_enst2type{$enst} = $type_label.":".$priority_n; - $priority_n++; - } - } - close(F); -} diff --git a/bin/perl/split_bam_to_subfiles_SEorPE.pl b/bin/perl/split_bam_to_subfiles_SEorPE.pl deleted file mode 100755 index df0707e..0000000 --- a/bin/perl/split_bam_to_subfiles_SEorPE.pl +++ /dev/null @@ -1,138 +0,0 @@ -#!/usr/bin/env perl - -use warnings; -use strict; - -my $sam_fi = $ARGV[0]; -my @sam_fi_split = split(/\//,$sam_fi); -my $sam_fi_short = $sam_fi_split[$#sam_fi_split]; - -my $se_or_pe_flag = $ARGV[1]; -if ($se_or_pe_flag eq "PE"){ - print STDERR "splitting in paired-end mode\n"; -} elsif ($se_or_pe_flag eq "SE") { - print STDERR "splitting in single-end mode\n"; -} else { - print STDERR "fatal error - SE or PE not defined\n"; - exit; -} - -my %filehandles; -for my $b1 ("A","C","G","T","N") { - for my $b2 ("A","C","G","T","N") { - - my $outfi = $b1.$b2.".".$sam_fi_short.".tmp"; - open(my $fh, '>', $outfi); - $filehandles{$b1.$b2} = $fh; - } -} - - - -if ($sam_fi =~ /\.sam$/) { - open(F,$sam_fi); -} elsif ($sam_fi =~ /\.bam$/) { - open(F,"samtools view -h $sam_fi |") || die "no $sam_fi\n"; -} else { - print STDERR "weird - $sam_fi not either sam or bam file format - exit\n"; - exit; -} - -while () { - my $r1 = $_; - - next if ($r1 =~ /^\@/); - - chomp($r1); - my @tmp_r1 = split(/\t/,$r1); - my ($r1name,$r1bc) = split(/\s+/,$tmp_r1[0]); - my $r1sam_flag = $tmp_r1[1]; - - my $r2 = if ($se_or_pe_flag eq "PE"); - chomp($r2) if ($se_or_pe_flag eq "PE"); - my @tmp_r2 = split(/\t/,$r2) if ($se_or_pe_flag eq "PE"); - my ($r2name,$r2bc) = split(/\s+/,$tmp_r2[0]) if ($se_or_pe_flag eq "PE"); - my $r2sam_flag = $tmp_r2[1] if ($se_or_pe_flag eq "PE"); - - if ($se_or_pe_flag eq "PE") { - unless ($r1name eq $r2name) { - print STDERR "paired end mismatch error: $sam_fi r1 $tmp_r1[0] r2 $tmp_r2[0]\n"; - } - unless ($r1sam_flag) { - print STDERR "error $r1 $r2\n"; - } - } - - - - if ($se_or_pe_flag eq "PE") { - next if ($r1sam_flag == 77 || $r1sam_flag == 141); - } elsif ($se_or_pe_flag eq "SE") { - next if ($r1sam_flag == 4); - } - my $frag_strand; -### This section is for only properly paired reads - if ($se_or_pe_flag eq "PE"){ - if ($r1sam_flag == 99 || $r1sam_flag == 355) { - $frag_strand = "-"; - } elsif ($r1sam_flag == 83 || $r1sam_flag == 339) { - $frag_strand = "+"; - } elsif ($r1sam_flag == 147 || $r1sam_flag == 403) { - $frag_strand = "-"; - @tmp_r1 = split(/\t/,$r2); - @tmp_r2 = split(/\t/,$r1); - } elsif ($r1sam_flag == 163 || $r1sam_flag == 419) { - $frag_strand = "+"; - @tmp_r1 = split(/\t/,$r2); - @tmp_r2 = split(/\t/,$r1); - } else { - next; - print STDERR "R1 strand error $r1sam_flag\n"; - } - my @read_name = split(/\:/,$tmp_r1[0]); - my $randommer = $read_name[0]; - my $first2rand = substr($randommer,0,2); - print { $filehandles{$first2rand} } "".$r1."\n".$r2."\n"; - - - - } elsif ($se_or_pe_flag eq "SE") { - if ($r1sam_flag == 16 || $r1sam_flag == 272) { - $frag_strand = "-"; - } elsif ($r1sam_flag eq "0" || $r1sam_flag == 256) { - $frag_strand = "+"; - } else { - next; - print STDERR "R1 strand error $r1sam_flag\n"; - } - - my @read_name = split(/\_/,$tmp_r1[0]); - my $randommer = pop(@read_name); - my $first2rand = substr($randommer,0,2); - print { $filehandles{$first2rand} } "".$r1."\n"; - - } -### - - # 77 = R1, unmapped - # 141 = R2, unmapped - # 99 = R1, mapped, fwd strand --- frag on rev strand -> 355 = not primary - # 147 = R2, mapped, rev strand -- frag on rev strand -> 403 = not primary - - # 101 = R1 unmapped, R2 mapped rev strand -- frag on rev strand - # 73 = R1, mapped, fwd strand --- frag on rev strand - # 153 = R2 mapped (R1 unmapped), rev strand -- frag on rev strand - # 133 = R2 unmapped, R1 mapped fwd strand -- frag on rev strand - - # 83 = R1, mapped, rev strand --- frag on fwd strand -> 339 = not primary - # 163 = R2, mapped, fwd strand -- frag on fwd strand -> 419 = not primary - - # 69 = R1 unmapped, R2 mapped fwd strand -- frag on fwd strand - # 89 = R1 mapped rev strand, R2 unmapped -- frag on fwd strand - # 137 = R2 mapped (R1 unmapped), fwd strand -- frag on fwd strand - # 165 = R2 unmapped, R1 rev strand -- frag on fwd strand - - - -} -close(F); diff --git a/bin/python/merge_multiple_parsed_files.simplified_20191022.py b/bin/python/merge_multiple_parsed_files.simplified_20191022.py new file mode 100755 index 0000000..42677a8 --- /dev/null +++ b/bin/python/merge_multiple_parsed_files.simplified_20191022.py @@ -0,0 +1,127 @@ +#!/usr/bin/env python3 +from __future__ import annotations + +import argparse +import fcntl +from collections import defaultdict +from pathlib import Path +import re +import sys + + +READINFO_RE = re.compile( + r'All\sreads\:\t(\d+)\tPCR\sduplicates\sremoved\:\t(\d+)\tUsable\sRemaining\:\t(\d+)\tUsable\sfrom\sgenomic\smapping\:\t(\d+)\tUsable\sfrom\sfamily\smapping\:\t(\d+)$' +) + + +def main() -> int: + ap = argparse.ArgumentParser(description='Merge parsed files') + ap.add_argument('output_fi') + ap.add_argument('files', nargs='+') + args = ap.parse_args() + + output_fi = Path(args.output_fi) + files = [Path(p) for p in args.files] + + split_fi1 = list(files[0].parts) + short_fi1 = split_fi1[-1] + working_dir = Path(*split_fi1[:-1]) if len(split_fi1) > 1 else Path('.') + failed_jobs_list = working_dir / 'failed_jobs_list.txt' + + print(f'SHORTFI:{short_fi1}', end='') + print(f'WORKDIR:{working_dir}', end='') + + lock_fi = Path(str(failed_jobs_list) + '.lck') + lock_fi.parent.mkdir(parents=True, exist_ok=True) + with lock_fi.open('w') as s: + fcntl.flock(s.fileno(), fcntl.LOCK_EX) + with failed_jobs_list.open('a'): + pass + + read_sums = { + 'all': 0, + 'usable': 0, + 'genomic': 0, + 'repfamily': 0, + 'total': defaultdict(int), + 'element': defaultdict(lambda: {'ensg_all': None, 'ensg_primary': None, 'readnum': 0}), + } + + for parsed_file in files: + parse_file(parsed_file, read_sums) + + output_fi.parent.mkdir(parents=True, exist_ok=True) + with output_fi.open('w', encoding='utf-8') as out: + usable = read_sums['usable'] + all_reads = read_sums['all'] + genomic = read_sums['genomic'] + repfamily = read_sums['repfamily'] + + out.write(f'#READINFO\tAllReads\t{all_reads}\n') + out.write(f'#READINFO\tUsableReads\t{usable}\t{usable / all_reads if all_reads else 0}\n') + out.write(f'#READINFO\tGenomicReads\t{genomic}\t{genomic / usable if usable else 0}\n') + out.write(f'#READINFO\tRepFamilyReads\t{repfamily}\t{repfamily / usable if usable else 0}\n') + + for element, readnum in sorted(read_sums['total'].items(), key=lambda kv: kv[1], reverse=True): + out.write(f'TOTAL\t{element}\t{readnum}\t{(readnum / usable) if usable else 0}\n') + + items = sorted( + read_sums['element'].items(), + key=lambda kv: kv[1]['readnum'], + reverse=True, + ) + for element, info in items: + out.write( + f"ELEMENT\t{info['ensg_primary']}\t{info['readnum']}\t" + f"{(info['readnum'] / usable) if usable else 0}\t{element}\t{info['ensg_all']}\n" + ) + + return 0 + + +def parse_file(parsed_fi: Path, read_sums: dict) -> None: + with parsed_fi.open('r', encoding='utf-8', errors='replace') as fh: + for raw in fh: + line = raw.rstrip('\n') + tmp = line.split('\t') + if not tmp: + continue + if tmp[0] == '#READINFO': + m = READINFO_RE.search(line) + if m: + read_sums['all'] += int(m.group(1)) + elif len(tmp) > 2 and tmp[1] == 'UsableReads': + read_sums['usable'] += int(float(tmp[2])) + elif len(tmp) > 2 and tmp[1] == 'GenomicReads': + read_sums['genomic'] += int(float(tmp[2])) + elif len(tmp) > 2 and tmp[1] == 'RepFamilyReads': + read_sums['repfamily'] += int(float(tmp[2])) + else: + print(f"couldn't parse readinfo line {line}", file=sys.stderr) + elif tmp[0] == '#READINFO2': + continue + elif tmp[0] == 'TOTAL': + _, element, readnum, *_ = tmp + read_sums['total'][element] += int(float(readnum)) + else: + if len(tmp) < 6: + continue + _, ensg_primary, readnum, _, enst_all, ensg_all = tmp[:6] + existing = read_sums['element'][enst_all] + if existing['ensg_all'] is not None and existing['ensg_all'] != ensg_all: + print( + f"error - ensg_all mismatch {ensg_all} {existing['ensg_all']}", + file=sys.stderr, + ) + if existing['ensg_primary'] is not None and existing['ensg_primary'] != ensg_primary: + print( + f"error - ensg_primary mismatch {ensg_primary} {existing['ensg_primary']}", + file=sys.stderr, + ) + existing['ensg_all'] = ensg_all + existing['ensg_primary'] = ensg_primary + existing['readnum'] += int(float(readnum)) + + +if __name__ == '__main__': + raise SystemExit(main()) diff --git a/bin/python/parse_bowtie2_output_realtime_includemultifamily_SE.py b/bin/python/parse_bowtie2_output_realtime_includemultifamily_SE.py new file mode 100755 index 0000000..920fd09 --- /dev/null +++ b/bin/python/parse_bowtie2_output_realtime_includemultifamily_SE.py @@ -0,0 +1,308 @@ +#!/usr/bin/env python3 +from __future__ import annotations + +import argparse +import subprocess +import sys +from collections import defaultdict +from dataclasses import dataclass, field +from pathlib import Path +from typing import Dict, List + + +RRNA_EXTRA_HASH = { + 'RNA28S': 'RNA45S', + 'RNA18S': 'RNA45S', + 'RNA5-8S': 'RNA45S', + 'antisense_RNA28S': 'antisense_RNA45S', + 'antisense_RNA18S': 'antisense_RNA45S', + 'antisense_RNA5-8S': 'antisense_RNA45S', +} +RRNA_EXTRA_HASH_REV = { + 'RNA45S': {'RNA28S', 'RNA18S', 'RNA5-8S'}, + 'antisense_RNA45S': {'antisense_RNA28S', 'antisense_RNA18S', 'antisense_RNA5-8S'}, +} + + +@dataclass +class ReadRecord: + r1: Dict[str, str] = field(default_factory=dict) + flags: Dict[str, int] = field(default_factory=dict) + quality: int = 0 + mult_ensts: Dict[str, List[str]] = field(default_factory=lambda: defaultdict(list)) + enst: Dict[str, str] = field(default_factory=dict) + master_enst: Dict[str, str] = field(default_factory=dict) + + +def read_in_filelists(filelist_file: Path) -> tuple[dict[str, str], dict[str, str]]: + enst2gene: dict[str, str] = {} + convert_enst2type: dict[str, str] = {} + priority_n = 0 + with filelist_file.open('r', encoding='utf-8', errors='replace') as fh: + for raw in fh: + line = raw.rstrip('\n') + if not line: + continue + parts = line.split('\t') + if len(parts) < 4: + continue + allenst, _allensg, gid, type_label = parts[:4] + type_label = type_label.rstrip('_') + if not allenst: + print(f'error missing enst {line} {filelist_file}', file=sys.stderr) + continue + for enst in allenst.split('|'): + enst2gene[enst] = gid + convert_enst2type[enst] = f'{type_label}:{priority_n}' + priority_n += 1 + return enst2gene, convert_enst2type + + +def parse_as_score(tags: list[str]) -> int: + for tag in tags: + if tag.startswith('AS:i:'): + try: + return int(tag[5:]) + except ValueError: + pass + return -10**9 + + +def print_output(read_hash: dict[str, ReadRecord], samout, multimapping_hash: dict[str, int]) -> None: + for read_name, rec in read_hash.items(): + ensttype_array = sorted(rec.flags.keys()) + if not ensttype_array: + continue + ensttype = ensttype_array[0] + masterenst_array = [rec.master_enst[t] for t in ensttype_array if t in rec.master_enst] + ensttype_join = '|'.join(ensttype_array) + masterenst_join = '|'.join(masterenst_array) + + if len(rec.flags) == 1: + r1_cols = rec.r1[ensttype].split('\t') + r1_cols[2] = f'{ensttype_join}||{masterenst_join}' + r1_line = '\t'.join(r1_cols) + zz = '|'.join(rec.mult_ensts.get(ensttype, [])) + samout.write(f'{r1_line}\tZZ:Z:{zz}\n') + else: + all_mult_ensts: list[str] = [] + for key in ensttype_array: + all_mult_ensts.append('|'.join(rec.mult_ensts.get(key, []))) + final_mult_ensts = '|'.join(all_mult_ensts) + if ensttype not in rec.r1: + print(f"weird error - {read_name} {ensttype} readhash doesn't exist ? {rec.flags.get(ensttype)}", file=sys.stderr) + continue + r1_cols = rec.r1[ensttype].split('\t') + r1_cols[2] = f'{ensttype_join}||{masterenst_join}' + r1_line = '\t'.join(r1_cols) + samout.write(f'{r1_line}\tZZ:Z:{final_mult_ensts}\n') + multimapping_type = '|'.join(rec.flags.keys()) + multimapping_hash[multimapping_type] = multimapping_hash.get(multimapping_type, 0) + 1 + + +def stream_sam_lines(args: argparse.Namespace): + if args.sam_input: + with Path(args.sam_input).open('r', encoding='utf-8', errors='replace') as fh: + for line in fh: + yield line + return + + bowtie_out = f'{args.output}.bowtieout' + cmd = [ + 'bowtie2', + '-q', + '--sensitive', + '-a', + '-p', + str(args.threads), + '--no-mixed', + '--reorder', + '-x', + args.bowtie_db, + '-U', + args.fastq_file1, + ] + print('command ' + ' '.join(cmd) + f' 2> {bowtie_out}', file=sys.stderr) + with Path(bowtie_out).open('w', encoding='utf-8') as err: + proc = subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=err, text=True) + assert proc.stdout is not None + for line in proc.stdout: + yield line + ret = proc.wait() + if ret != 0: + raise RuntimeError(f'bowtie2 failed with exit code {ret}') + + +def main() -> int: + ap = argparse.ArgumentParser(description='Python port of parse_bowtie2_output_realtime_includemultifamily_SE.pl') + ap.add_argument('fastq_file1') + ap.add_argument('bowtie_db') + ap.add_argument('output') + ap.add_argument('filelist_file') + ap.add_argument('--sam-input', default=None, help='Optional SAM input stream file for testing/debugging') + ap.add_argument('--threads', type=int, default=3) + args = ap.parse_args() + + _enst2gene, convert_enst2type = read_in_filelists(Path(args.filelist_file)) + + output_path = Path(args.output) + multimapping_out = Path(args.output + '.multimapping_deleted') + done_file = Path(args.output + '.done') + + read_hash: dict[str, ReadRecord] = {} + multimapping_hash: dict[str, int] = {} + print_batch = 10000 + read_counter = 0 + prev_r1name = '' + + with output_path.open('w', encoding='utf-8') as samout: + for raw in stream_sam_lines(args): + r1 = raw.rstrip('\n') + if r1.startswith('@'): + samout.write(r1 + '\n') + continue + + cols = r1.split('\t') + if len(cols) < 12: + continue + + read_name_parts = cols[0].split('_') + if len(read_name_parts) > 1: + r1name = '_'.join(read_name_parts[:-1]) + else: + r1name = cols[0] + + try: + r1sam_flag = int(cols[1]) + except ValueError: + continue + if r1sam_flag == 4: + continue + + if r1sam_flag in (16, 272): + frag_strand = '-' + elif r1sam_flag in (0, 256): + frag_strand = '+' + else: + continue + + paired_mismatch_score = parse_as_score(cols[11:]) + + mapped_enst = cols[2] + mapped_enst_full = cols[2] + if mapped_enst.endswith('_spliced'): + mapped_enst = mapped_enst[: -len('_spliced')] + if mapped_enst.endswith('_withgenomeflank'): + mapped_enst = mapped_enst[: -len('_withgenomeflank')] + + if mapped_enst not in convert_enst2type: + print(f'enst2type is missing for {mapped_enst} {r1}', file=sys.stderr) + continue + ensttype, enstpriority_s = convert_enst2type[mapped_enst].split(':', 1) + enstpriority = int(enstpriority_s) + + if frag_strand == '-': + ensttype = 'antisense_' + ensttype + mapped_enst_full = 'antisense_' + mapped_enst_full + + if r1name != prev_r1name: + read_counter += 1 + if read_counter > print_batch: + print_output(read_hash, samout, multimapping_hash) + read_hash = {} + read_counter = 0 + prev_r1name = r1name + + if r1name not in read_hash: + rec = ReadRecord() + rec.r1[ensttype] = r1 + rec.flags[ensttype] = enstpriority + rec.quality = paired_mismatch_score + rec.mult_ensts[ensttype].append(mapped_enst_full) + rec.enst[mapped_enst] = mapped_enst_full + rec.master_enst[ensttype] = mapped_enst_full + read_hash[r1name] = rec + continue + + rec = read_hash[r1name] + if paired_mismatch_score < rec.quality: + continue + if paired_mismatch_score > rec.quality: + newrec = ReadRecord() + newrec.r1[ensttype] = r1 + newrec.flags[ensttype] = enstpriority + newrec.quality = paired_mismatch_score + newrec.mult_ensts[ensttype].append(mapped_enst_full) + newrec.enst[mapped_enst] = mapped_enst_full + newrec.master_enst[ensttype] = mapped_enst_full + read_hash[r1name] = newrec + continue + + if ensttype in rec.flags: + if mapped_enst in rec.enst: + old_full = rec.enst[mapped_enst] + if old_full + '_withgenomeflank' == mapped_enst_full or old_full + '_spliced' == mapped_enst_full: + pass + elif old_full == mapped_enst_full + '_withgenomeflank' or old_full == mapped_enst_full + '_spliced': + rec.r1[ensttype] = r1 + rec.flags[ensttype] = enstpriority + for i, v in enumerate(rec.mult_ensts.get(ensttype, [])): + if v == old_full: + rec.mult_ensts[ensttype][i] = mapped_enst_full + rec.enst[mapped_enst] = mapped_enst_full + rec.master_enst[ensttype] = mapped_enst_full + else: + for i, v in enumerate(rec.mult_ensts.get(ensttype, [])): + if v == old_full: + rec.mult_ensts[ensttype][i] = mapped_enst_full + '_DOUBLEMAP' + rec.master_enst[ensttype] = mapped_enst_full + '_DOUBLEMAP' + elif enstpriority < rec.flags[ensttype]: + rec.r1[ensttype] = r1 + rec.flags[ensttype] = enstpriority + rec.mult_ensts[ensttype].insert(0, mapped_enst_full) + rec.enst[mapped_enst] = mapped_enst_full + rec.master_enst[ensttype] = mapped_enst_full + else: + rec.mult_ensts[ensttype].append(mapped_enst_full) + elif ensttype in RRNA_EXTRA_HASH and RRNA_EXTRA_HASH[ensttype] in rec.r1: + old_rrna = RRNA_EXTRA_HASH[ensttype] + rec.r1.pop(old_rrna, None) + rec.flags.pop(old_rrna, None) + rec.master_enst.pop(old_rrna, None) + rec.mult_ensts.pop(old_rrna, None) + rec.enst.pop('NR_046235.1', None) + rec.r1[ensttype] = r1 + rec.flags[ensttype] = enstpriority + rec.quality = paired_mismatch_score + rec.enst[mapped_enst] = mapped_enst_full + rec.master_enst[ensttype] = mapped_enst_full + rec.mult_ensts[ensttype].append(mapped_enst_full) + elif ensttype in RRNA_EXTRA_HASH_REV: + rna_flag = any(el in rec.r1 for el in RRNA_EXTRA_HASH_REV[ensttype]) + if not rna_flag: + rec.r1[ensttype] = r1 + rec.flags[ensttype] = enstpriority + rec.quality = paired_mismatch_score + rec.enst[mapped_enst] = mapped_enst_full + rec.mult_ensts[ensttype].append(mapped_enst_full) + rec.master_enst[ensttype] = mapped_enst_full + else: + rec.r1[ensttype] = r1 + rec.flags[ensttype] = enstpriority + rec.quality = paired_mismatch_score + rec.enst[mapped_enst] = mapped_enst_full + rec.mult_ensts[ensttype].append(mapped_enst_full) + rec.master_enst[ensttype] = mapped_enst_full + + print_output(read_hash, samout, multimapping_hash) + + with multimapping_out.open('w', encoding='utf-8') as mm: + for key, count in multimapping_hash.items(): + mm.write(f'{key}\t{count}\n') + + done_file.write_text('jobs done\n', encoding='utf-8') + return 0 + + +if __name__ == '__main__': + raise SystemExit(main()) diff --git a/bin/python/split_bam_to_subfiles_SEorPE.py b/bin/python/split_bam_to_subfiles_SEorPE.py new file mode 100755 index 0000000..0b4c237 --- /dev/null +++ b/bin/python/split_bam_to_subfiles_SEorPE.py @@ -0,0 +1,120 @@ +#!/usr/bin/env python3 +from __future__ import annotations + +import argparse +import subprocess +import sys +from pathlib import Path + + +PREFIXES = [a + b for a in 'ACGTN' for b in 'ACGTN'] + + +def sam_stream(path: Path): + if path.suffix == '.sam': + with path.open('r', encoding='utf-8', errors='replace') as fh: + for line in fh: + yield line + elif path.suffix == '.bam': + proc = subprocess.Popen( + ['samtools', 'view', '-h', str(path)], + stdout=subprocess.PIPE, + stderr=subprocess.PIPE, + text=True, + ) + assert proc.stdout is not None + for line in proc.stdout: + yield line + stderr = proc.stderr.read() if proc.stderr else '' + ret = proc.wait() + if ret != 0: + raise RuntimeError(f'samtools view failed ({ret}): {stderr}') + else: + raise ValueError(f'weird - {path} not either sam or bam file format - exit') + + +def main() -> int: + ap = argparse.ArgumentParser(description='Split SAM/BAM into 25 UMI-prefix tmp files') + ap.add_argument('sam_fi') + ap.add_argument('se_or_pe_flag', choices=['SE', 'PE']) + args = ap.parse_args() + + sam_path = Path(args.sam_fi) + sam_short = sam_path.name + mode = args.se_or_pe_flag + + if mode == 'PE': + print('splitting in paired-end mode', file=sys.stderr) + else: + print('splitting in single-end mode', file=sys.stderr) + + filehandles = { + p: (Path(f'{p}.{sam_short}.tmp')).open('w', encoding='utf-8') + for p in PREFIXES + } + + try: + stream = sam_stream(sam_path) + it = iter(stream) + for r1 in it: + if r1.startswith('@'): + continue + r1 = r1.rstrip('\n') + tmp_r1 = r1.split('\t') + r1_name = tmp_r1[0].split()[0] + r1_flag = int(tmp_r1[1]) + + r2 = None + tmp_r2: list[str] | None = None + if mode == 'PE': + try: + r2 = next(it).rstrip('\n') + except StopIteration: + break + tmp_r2 = r2.split('\t') + r2_name = tmp_r2[0].split()[0] + if r1_name != r2_name: + print( + f'paired end mismatch error: {sam_path} r1 {tmp_r1[0]} r2 {tmp_r2[0]}', + file=sys.stderr, + ) + if not r1_flag: + print(f'error {r1} {r2}', file=sys.stderr) + + if mode == 'PE': + if r1_flag in (77, 141): + continue + else: + if r1_flag == 4: + continue + + if mode == 'PE': + if r1_flag in (99, 355, 83, 339): + pass + elif r1_flag in (147, 403, 163, 419): + assert r2 is not None and tmp_r2 is not None + tmp_r1, tmp_r2 = r2.split('\t'), r1.split('\t') + else: + continue + + randommer = tmp_r1[0].split(':')[0] + first2 = randommer[:2] + if first2 in filehandles and r2 is not None: + filehandles[first2].write(f'{r1}\n{r2}\n') + else: + if r1_flag not in (16, 272, 0, 256): + continue + read_name_parts = tmp_r1[0].split('_') + randommer = read_name_parts[-1] + first2 = randommer[:2] + if first2 in filehandles: + filehandles[first2].write(f'{r1}\n') + finally: + for fh in filehandles.values(): + fh.close() + + return 0 + + +if __name__ == '__main__': + raise SystemExit(main()) diff --git a/changelog.md b/changelog.md index 28b7bd6..050d5a7 100644 --- a/changelog.md +++ b/changelog.md @@ -48,3 +48,6 @@ The format is based on [Keep a Changelog](http://keepachangelog.com/en/1.0.0/) - Generates repeat-mapped rmDup and preRmDup SAM-like files mapping reads to repeat element families. - README detailing methods and output explanations + +## codex/python-conversion +- removed legacy NOTES file; branch instructions now centralized in README.md diff --git a/config/config.yaml b/config/config.yaml new file mode 100644 index 0000000..004a4c0 --- /dev/null +++ b/config/config.yaml @@ -0,0 +1,50 @@ +# CWL-compatible input names for SE pipeline wiring +# Fill these with real dataset paths for production runs. +dataset: seCLIP_example +barcode1r1FastqGz: "" +barcode1rmRepBam: "" +barcode1Inputr1FastqGz: "" +barcode1InputrmRepBam: "" + +bowtie2_db: "" +bowtie2_prefix: "" +fileListFile1: "" +gencodeGTF: "" +gencodeTableBrowser: "" +repMaskBEDFile: "" + +prefixes: + - AA + - AC + - AG + - AT + - AN + - CA + - CC + - CG + - CT + - CN + - GA + - GC + - GG + - GT + - GN + - TA + - TC + - TG + - TT + - TN + - NA + - NC + - NG + - NT + - NN +se_or_pe: SE + +# Mini validation fixture configuration (used by currently implemented rules) +mini_validation: + source_sam_gz: tests/fixtures/mini/source/ip.preRmDup.sam.mini.gz + split_manifest: tests/fixtures/mini/expected/split.expected.manifest.tsv + merge_input_1: tests/fixtures/mini/source/merge_input_1.parsed_v2.txt + merge_input_2: tests/fixtures/mini/source/merge_input_2.parsed_v2.txt + merge_expected: tests/fixtures/mini/expected/merged.expected.parsed diff --git a/cwl/calculate_fold_change_from_parsed_files.cwl b/cwl/calculate_fold_change_from_parsed_files.cwl deleted file mode 100755 index 9ed2203..0000000 --- a/cwl/calculate_fold_change_from_parsed_files.cwl +++ /dev/null @@ -1,88 +0,0 @@ -#!/usr/bin/env cwltool - -cwlVersion: v1.0 - -class: CommandLineTool - -requirements: - - class: InlineJavascriptRequirement - -hints: - - class: DockerRequirement - dockerPull: brianyee/repetitive_element_mapping:1.0.0 - -baseCommand: [calculate_fold_change_from_parsed_files.py] - -inputs: - - ip_parsed_file: - type: File - inputBinding: - prefix: --ip_parsed - position: 1 - - input_parsed_file: - type: File - inputBinding: - prefix: --input_parsed - position: 2 - - out_file_nopipes: - type: string - inputBinding: - position: 3 - prefix: --out_file_nopipes - valueFrom: | - ${ - if (inputs.out_file_nopipes == "") { - return inputs.ip_parsed_file.nameroot + ".nopipes.tsv"; - } - else { - return inputs.out_file_nopipes; - } - } - default: "" - - out_file_withpipes: - type: string - inputBinding: - position: 4 - prefix: --out_file_withpipes - valueFrom: | - ${ - if (inputs.out_file_nopipes == "") { - return inputs.ip_parsed_file.nameroot + ".withpipes.tsv"; - } - else { - return inputs.out_file_withpipes; - } - } - default: "" - -outputs: - - out_file_nopipes_file: - type: File - outputBinding: - glob: | - ${ - if (inputs.out_file_nopipes == "") { - return inputs.ip_parsed_file.nameroot + ".nopipes.tsv"; - } - else { - return inputs.out_file_nopipes; - } - } - - out_file_withpipes_file: - type: File - outputBinding: - glob: | - ${ - if (inputs.out_file_withpipes == "") { - return inputs.ip_parsed_file.nameroot + ".withpipes.tsv"; - } - else { - return inputs.out_file_withpipes; - } - } \ No newline at end of file diff --git a/cwl/combine.cwl b/cwl/combine.cwl deleted file mode 100755 index cb9a3f2..0000000 --- a/cwl/combine.cwl +++ /dev/null @@ -1,37 +0,0 @@ -#!/usr/bin/env cwltool - -cwlVersion: v1.0 - -class: CommandLineTool - -requirements: - - class: InlineJavascriptRequirement - - class: InitialWorkDirRequirement - listing: $(inputs.files) - -hints: - - class: DockerRequirement - dockerPull: brianyee/repetitive_element_mapping:1.0.0 - -baseCommand: [merge_multiple_parsed_files.simplified_20191022.pl] - -inputs: - - outputFile: - type: string - inputBinding: - position: 1 - - files: - type: File[] - inputBinding: - position: 2 - -outputs: - - output: - type: File - outputBinding: - glob: $(inputs.outputFile) - - diff --git a/cwl/concatenate.cwl b/cwl/concatenate.cwl deleted file mode 100755 index 67b1b13..0000000 --- a/cwl/concatenate.cwl +++ /dev/null @@ -1,26 +0,0 @@ -#!/usr/bin/env cwltool - -cwlVersion: v1.0 - -class: CommandLineTool - -baseCommand: [cat] - -inputs: - - files: - type: File[] - inputBinding: - position: 1 - - concatenated_output: - type: string - -stdout: $(inputs.concatenated_output) - -outputs: - - concatenated: - type: File - outputBinding: - glob: $(inputs.concatenated_output) diff --git a/cwl/deduplicate.cwl b/cwl/deduplicate.cwl deleted file mode 100755 index c07e8e6..0000000 --- a/cwl/deduplicate.cwl +++ /dev/null @@ -1,90 +0,0 @@ -#!/usr/bin/env cwltool - -cwlVersion: v1.0 - -class: CommandLineTool - -requirements: - - class: ResourceRequirement - ramMin: 32000 - -hints: - - class: DockerRequirement - dockerPull: brianyee/repetitive_element_mapping:1.0.0 - -baseCommand: [duplicate_removal_inline_paired.count_region_other_reads_masksnRNAs_andreparse_SEandPE_20201210_simple.pl] - -inputs: - - repFamilySam: - type: File - inputBinding: - position: 1 - label: "input" - doc: "input sam file" - - rmRepSam: - type: File - inputBinding: - position: 2 - label: "rmRep.sam alignment file" - doc: "rmRep sam alignment file" - - se_or_pe: - type: string - inputBinding: - position: 3 - doc: "Either PE or SE (uppercase required)" - - gencodeGTF: - type: File - inputBinding: - position: 4 - label: "gencode GTF File" - doc: "gencode GTF File" - - gencodeTableBrowser: - type: File - inputBinding: - position: 5 - label: "gencode GTF file in UCSC table browser format" - doc: "gencode GTF file in UCSC table browser format" - - repMaskBedFile: - type: File - inputBinding: - position: 6 - label: "repeatmasker bedfile" - doc: "bedfile of repeat regions" - - fileList1: - type: File - inputBinding: - position: 7 - doc: "tsv with 5 fields: ENST/ENSG/name/chrom/genelist.file" - -outputs: - - deduplicatedRmDupSam: - type: File - outputBinding: - glob: "*combined_w_uniquemap.rmDup.sam" - label: "combined unique mapping + rep element mapping sam" - - deduplicatedPreRmDupSam: - type: File - outputBinding: - glob: "*combined_w_uniquemap.prermDup.sam" - label: "combined unique mapping + rep element mapping sam before rmdup" - - parsedFile: - type: File - outputBinding: - glob: "*.parsed_v2.20201210.txt" - label: "combined unique mapping + rep element mapping sam" - - doneFile: - type: File - outputBinding: - glob: "*.done" - label: "done file" \ No newline at end of file diff --git a/cwl/getpair.cwl b/cwl/getpair.cwl deleted file mode 100755 index e68440c..0000000 --- a/cwl/getpair.cwl +++ /dev/null @@ -1,50 +0,0 @@ -#!/usr/bin/env cwltool - -cwlVersion: v1.0 - -class: ExpressionTool - -requirements: - - class: InlineJavascriptRequirement - - class: ResourceRequirement - ramMin: 1000 - -inputs: - - prefix: - type: string - doc: 'First two letters of barcode ID (AA, AC, ..., NN)' - rep_s: - type: File[] - doc: 'List of repetitive element-mapped sam-like files, split by prefix' - rmrep_s: - type: File[] - doc: 'List of remove-rep BAM files, split by prefix' - -outputs: - - prefixrep: - type: File - prefixrmrep: - type: File - -expression: | - ${ - var prefix = inputs.prefix; - var rep_s = inputs.rep_s; - var rmrep_s = inputs.rmrep_s; - - var prefixrep = ''; - var prefixrmrep = ''; - - for (var i = 0; i < rep_s.length; i++) { - if (rep_s[i].basename.indexOf(prefix) == 0) { - prefixrep = rep_s[i]; - } - if (rmrep_s[i].basename.indexOf(prefix) == 0) { - prefixrmrep = rmrep_s[i]; - } - } - return {'prefixrep': prefixrep, 'prefixrmrep': prefixrmrep} - } - diff --git a/cwl/gzip.cwl b/cwl/gzip.cwl deleted file mode 100755 index 4748cd5..0000000 --- a/cwl/gzip.cwl +++ /dev/null @@ -1,30 +0,0 @@ -#!/usr/bin/env cwltool - -cwlVersion: v1.0 - -class: CommandLineTool - -baseCommand: [gzip] - -inputs: - - stdout: - type: boolean - inputBinding: - position: 1 - prefix: -c - default: true - - input: - type: File - inputBinding: - position: 2 - -stdout: $(inputs.input.basename).gz - -outputs: - - gzipped: - type: File - outputBinding: - glob: $(inputs.input.basename).gz diff --git a/cwl/map_repetitive_elements_pe.cwl b/cwl/map_repetitive_elements_pe.cwl deleted file mode 100755 index 3391f22..0000000 --- a/cwl/map_repetitive_elements_pe.cwl +++ /dev/null @@ -1,80 +0,0 @@ -#!/usr/bin/env cwltool - -cwlVersion: v1.0 - -class: CommandLineTool - -requirements: - - class: InlineJavascriptRequirement - - class: ResourceRequirement - coresMin: 8 - ramMin: 16000 - -hints: - - class: DockerRequirement - dockerPull: brianyee/repetitive_element_mapping:1.0.0 - -# # wrapped perl script for parsing bowtie results inline -baseCommand: [parse_bowtie2_output_realtime_includemultifamily_PE.pl] - - -inputs: - - - #read1 trimmed fastq - read1: - type: File - inputBinding: - position: 1 - - # read2 trimmed fastq" - read2: - type: File - inputBinding: - position: 2 - - bowtie2_db: - type: Directory - - bowtie2_prefix: - type: string - - output_file: - default: "" - type: string - inputBinding: - position: 4 - valueFrom: | - ${ - if (inputs.output_file == "") { - return inputs.read1.nameroot + ".Rep.sam"; - } - else { - return inputs.output_file; - } - } - file_list_file: - type: File - inputBinding: - position: 5 - -arguments: - - valueFrom: $(inputs.bowtie2_db.path + "//" + inputs.bowtie2_prefix) - position: 3 - shellQuote: false - -outputs: - - rep_sam: - type: File - outputBinding: - glob: | - ${ - if (inputs.output_file == "") { - return inputs.read1.nameroot + ".Rep.sam"; - } - else { - return inputs.output_file; - } - } - diff --git a/cwl/map_repetitive_elements_se.cwl b/cwl/map_repetitive_elements_se.cwl deleted file mode 100755 index e386316..0000000 --- a/cwl/map_repetitive_elements_se.cwl +++ /dev/null @@ -1,73 +0,0 @@ -#!/usr/bin/env cwltool - -cwlVersion: v1.0 - -class: CommandLineTool - -requirements: - - class: InlineJavascriptRequirement - - class: ResourceRequirement - coresMin: 8 - ramMin: 16000 - -hints: - - class: DockerRequirement - dockerPull: brianyee/repetitive_element_mapping:1.0.0 - -# # wrapped perl script for parsing bowtie results inline -baseCommand: [parse_bowtie2_output_realtime_includemultifamily_SE.pl] - - -inputs: - - - #read1 trimmed fastq - read1: - type: File - inputBinding: - position: 1 - - bowtie2_db: - type: Directory - - bowtie2_prefix: - type: string - - output_file: - default: "" - type: string - inputBinding: - position: 3 - valueFrom: | - ${ - if (inputs.output_file == "") { - return inputs.read1.nameroot + ".Rep.sam"; - } - else { - return inputs.output_file; - } - } - file_list_file: - type: File - inputBinding: - position: 4 - -arguments: - - valueFrom: $(inputs.bowtie2_db.path + "//" + inputs.bowtie2_prefix) - position: 2 - shellQuote: false - -outputs: - - rep_sam: - type: File - outputBinding: - glob: | - ${ - if (inputs.output_file == "") { - return inputs.read1.nameroot + ".Rep.sam"; - } - else { - return inputs.output_file; - } - } diff --git a/cwl/splitbam.cwl b/cwl/splitbam.cwl deleted file mode 100755 index 445fc62..0000000 --- a/cwl/splitbam.cwl +++ /dev/null @@ -1,35 +0,0 @@ -#!/usr/bin/env cwltool - -cwlVersion: v1.0 - -class: CommandLineTool - -hints: - - class: DockerRequirement - dockerPull: brianyee/repetitive_element_mapping:1.0.0 - -baseCommand: [split_bam_to_subfiles_SEorPE.pl] - -requirements: - - class: InlineJavascriptRequirement - -inputs: - - sam_file: - type: File - inputBinding: - position: 1 - se_or_pe: - type: string - inputBinding: - position: 2 - doc: "Either PE or SE (uppercase required)" - -outputs: - - repsam_s: - type: - type: array - items: File - outputBinding: - glob: "*.tmp" diff --git a/cwl/wf_ecliprepmap_pe.cwl b/cwl/wf_ecliprepmap_pe.cwl deleted file mode 100755 index b248163..0000000 --- a/cwl/wf_ecliprepmap_pe.cwl +++ /dev/null @@ -1,278 +0,0 @@ -#!/usr/bin/env cwltool - -cwlVersion: v1.0 -class: Workflow - -requirements: - - class: SubworkflowFeatureRequirement - - class: MultipleInputFeatureRequirement - - class: InlineJavascriptRequirement - - class: StepInputExpressionRequirement - -inputs: - - barcode1r1FastqGz: - type: File - barcode1r2FastqGz: - type: File - barcode1rmRepBam: - type: File - - barcode2r1FastqGz: - type: File - barcode2r2FastqGz: - type: File - barcode2rmRepBam: - type: File - - barcode1Inputr1FastqGz: - type: File - barcode1Inputr2FastqGz: - type: File - barcode1InputrmRepBam: - type: File - - bowtie2_db: - type: Directory - bowtie2_prefix: - type: string - fileListFile1: - type: File - - gencodeGTF: - type: File - gencodeTableBrowser: - type: File - repMaskBEDFile: - type: File - - prefixes: - type: string[] - default: [ - "AA","AC","AG","AT","AN", - "CA","CC","CG","CT","CN", - "GA","GC","GG","GT","GN", - "TA","TC","TG","TT","TN", - "NA","NC","NG","NT","NN" - ] - se_or_pe: - type: string - default: PE - -outputs: - - - ### PRE RMDUPED SAM FILE FINAL OUTPUTS ### - - - output_ip_concatenated_pre_rmDup_sam_file: - type: File - outputSource: step_gzip_preRmDup/gzipped - output_input_concatenated_pre_rmDup_sam_file: - type: File - outputSource: step_ecliprepmap_input/output_concatenated_preRmDup_sam_file - - - ### RMDUPED SAM FILE FINAL OUTPUTS ### - - output_barcode1_concatenated_rmDup_sam_file: - type: File - outputSource: step_gzip_rmDup/gzipped - output_input_concatenated_rmDup_sam_file: - type: File - outputSource: step_ecliprepmap_input/output_concatenated_rmDup_sam_file - - - ### FINAL PARSED STATISTICS FILES ### - - output_ip_parsed: - type: File - outputSource: step_combine_parsed/output - output_input_parsed: - type: File - outputSource: step_ecliprepmap_input/output_combined_parsed_file - output_nopipes: - type: File - outputSource: step_calculate_fold_change_from_parsed_files/out_file_nopipes_file - output_withpipes: - type: File - outputSource: step_calculate_fold_change_from_parsed_files/out_file_withpipes_file - -steps: - -########################################################################### -# Repeat-map paired-end IP (2 barcodes) -########################################################################### - - step_ecliprepmap_barcode1: - run: wf_ecliprepmap_pe_1barcode.cwl - in: - dataset: - source: barcode1r1FastqGz - valueFrom: | - ${ - return self.nameroot + ".barcode1"; - } - r1FastqGz: barcode1r1FastqGz - r2FastqGz: barcode1r2FastqGz - rmRepBam: barcode1rmRepBam - bowtie2_db: bowtie2_db - bowtie2_prefix: bowtie2_prefix - fileListFile1: fileListFile1 - gencodeGTF: gencodeGTF - gencodeTableBrowser: gencodeTableBrowser - repMaskBEDFile: repMaskBEDFile - prefixes: prefixes - se_or_pe: se_or_pe - out: - - output_repeat_mapped_sam_file - - output_rmDup_sam_files - - output_pre_rmDup_sam_files - - output_concatenated_rmDup_sam_file - - output_concatenated_preRmDup_sam_file - - output_parsed_files - - output_combined_parsed_file - - step_ecliprepmap_barcode2: - run: wf_ecliprepmap_pe_1barcode.cwl - in: - dataset: - source: barcode2r1FastqGz - valueFrom: | - ${ - return self.nameroot + ".barcode2"; - } - r1FastqGz: barcode2r1FastqGz - r2FastqGz: barcode2r2FastqGz - rmRepBam: barcode2rmRepBam - bowtie2_db: bowtie2_db - bowtie2_prefix: bowtie2_prefix - fileListFile1: fileListFile1 - gencodeGTF: gencodeGTF - gencodeTableBrowser: gencodeTableBrowser - repMaskBEDFile: repMaskBEDFile - prefixes: prefixes - se_or_pe: se_or_pe - out: - - output_repeat_mapped_sam_file - - output_rmDup_sam_files - - output_pre_rmDup_sam_files - - output_concatenated_rmDup_sam_file - - output_concatenated_preRmDup_sam_file - - output_parsed_files - - output_combined_parsed_file - -########################################################################### -# Repeat-map input sample (1 sample) -########################################################################### - - step_ecliprepmap_input: - run: wf_ecliprepmap_pe_1barcode.cwl - in: - dataset: - source: barcode1Inputr1FastqGz - valueFrom: | - ${ - return self.nameroot + ".input"; - } - r1FastqGz: barcode1Inputr1FastqGz - r2FastqGz: barcode1Inputr2FastqGz - bowtie2_db: bowtie2_db - bowtie2_prefix: bowtie2_prefix - fileListFile1: fileListFile1 - rmRepBam: barcode1InputrmRepBam - gencodeGTF: gencodeGTF - gencodeTableBrowser: gencodeTableBrowser - repMaskBEDFile: repMaskBEDFile - prefixes: prefixes - se_or_pe: se_or_pe - out: - - output_repeat_mapped_sam_file - - output_rmDup_sam_files - - output_pre_rmDup_sam_files - - output_concatenated_rmDup_sam_file - - output_concatenated_preRmDup_sam_file - - output_parsed_files - - output_combined_parsed_file - -########################################################################### -# Combine parsed files and calculate fold change/entropy -########################################################################### - - step_combine_parsed: - doc: "concatenates all final statistics using custom perl script" - run: combine.cwl - in: - files: - source: [ - step_ecliprepmap_barcode1/output_parsed_files, - step_ecliprepmap_barcode2/output_parsed_files - ] - linkMerge: merge_flattened - outputFile: - source: barcode1r1FastqGz - valueFrom: | - ${ - return self.nameroot + ".combined.parsed"; - } - out: - - output - - step_calculate_fold_change_from_parsed_files: - run: calculate_fold_change_from_parsed_files.cwl - in: - ip_parsed_file: step_combine_parsed/output - input_parsed_file: step_ecliprepmap_input/output_combined_parsed_file - out: - - out_file_nopipes_file - - out_file_withpipes_file - - step_concatenate_rmDup: - doc: "concatenates all rmdup sam files using cat" - run: concatenate.cwl - in: - files: - source: - - step_ecliprepmap_barcode1/output_rmDup_sam_files - - step_ecliprepmap_barcode2/output_rmDup_sam_files - linkMerge: merge_flattened - concatenated_output: - source: barcode1r1FastqGz - valueFrom: | - ${ - return self.nameroot + ".rmDup.sam"; - } - out: - - concatenated - - step_gzip_rmDup: - run: gzip.cwl - in: - input: step_concatenate_rmDup/concatenated - out: - - gzipped - - step_concatenate_pre_rmDup: - doc: "concatenates all pre-rmduped sam files using cat" - run: concatenate.cwl - in: - files: - source: - - step_ecliprepmap_barcode1/output_pre_rmDup_sam_files - - step_ecliprepmap_barcode2/output_pre_rmDup_sam_files - linkMerge: merge_flattened - concatenated_output: - source: barcode1r1FastqGz - valueFrom: | - ${ - return self.nameroot + ".preRmDup.sam"; - } - out: - - concatenated - - step_gzip_preRmDup: - run: gzip.cwl - in: - input: step_concatenate_pre_rmDup/concatenated - out: - - gzipped diff --git a/cwl/wf_ecliprepmap_pe_1barcode.cwl b/cwl/wf_ecliprepmap_pe_1barcode.cwl deleted file mode 100755 index 0a67952..0000000 --- a/cwl/wf_ecliprepmap_pe_1barcode.cwl +++ /dev/null @@ -1,213 +0,0 @@ -#!/usr/bin/env cwltool - -cwlVersion: v1.0 -class: Workflow - -requirements: - - class: ScatterFeatureRequirement - - class: InlineJavascriptRequirement - - class: StepInputExpressionRequirement - -inputs: - - dataset: - type: string - - r1FastqGz: - doc: "Read1 after trimming inline barcodes" - type: File - r2FastqGz: - doc: "Read2 after trimming inline barcodes" - type: File - rmRepBam: - doc: "BAM file after removing repetitive elements (after 2nd STAR mapping)" - type: File - - bowtie2_db: - type: Directory - bowtie2_prefix: - type: string - - fileListFile1: - type: File - - gencodeGTF: - type: File - gencodeTableBrowser: - type: File - repMaskBEDFile: - type: File - - prefixes: - type: string[] - default: [ - "AA","AC","AG","AT","AN", - "CA","CC","CG","CT","CN", - "GA","GC","GG","GT","GN", - "TA","TC","TG","TT","TN", - "NA","NC","NG","NT","NN" - ] - - se_or_pe: - type: string - default: PE - -outputs: - - output_repeat_mapped_sam_file: - doc: "first mapped SAM-like file to bowtie2" - type: File - outputSource: step_map_repetitive_elements/rep_sam - - output_rmDup_sam_files: - doc: "duplicate removed SAM-like files, one for each barcode" - type: File[] - outputSource: step_deduplicate/deduplicatedRmDupSam - - output_pre_rmDup_sam_files: - doc: "pre-duplicate removed SAM-like files, one for each barcode" - type: File[] - outputSource: step_deduplicate/deduplicatedPreRmDupSam - - output_concatenated_rmDup_sam_file: - doc: "Final remove duplicate SAM-like file" - type: File - outputSource: step_gzip_rmDup/gzipped - - output_concatenated_preRmDup_sam_file: - doc: "Final duplicated SAM-like file" - type: File - outputSource: step_gzip_preRmDup/gzipped - - output_parsed_files: - doc: "Output file containing read stats for each UMI prefix" - type: File[] - outputSource: step_deduplicate/parsedFile - - output_combined_parsed_file: - doc: "Combined output file containing read stats for all UMI prefixes" - type: File - outputSource: step_combine_parsed/output - -steps: - - step_map_repetitive_elements: - run: map_repetitive_elements_pe.cwl - in: - read1: r1FastqGz - read2: r2FastqGz - bowtie2_db: bowtie2_db - bowtie2_prefix: bowtie2_prefix - file_list_file: fileListFile1 - out: - - rep_sam - - step_splitbam_repsam: - run: splitbam.cwl - in: - sam_file: step_map_repetitive_elements/rep_sam - se_or_pe: se_or_pe - out: - - repsam_s - - step_splitbam_rmrepbam: - run: splitbam.cwl - in: - sam_file: rmRepBam - se_or_pe: se_or_pe - out: - - repsam_s - - step_getpair: - doc: | - Given a prefix (AA, AC, ... NN), return the rep and rmrep pairs - belonging to each prefix. Each file contains reads whose first 2nt - of its umi matches the prefix. - run: getpair.cwl - in: - rep_s: step_splitbam_repsam/repsam_s - rmrep_s: step_splitbam_rmrepbam/repsam_s - prefix: prefixes - scatter: prefix - out: - - prefixrep - - prefixrmrep - - step_deduplicate: - doc: | - Takes the repeat-mapped sam-like file, and the remove-replicate - bam file, and uses UMI information to remove PCR duplicates. - run: deduplicate.cwl - in: - repFamilySam: step_getpair/prefixrep - rmRepSam: step_getpair/prefixrmrep - se_or_pe: se_or_pe - gencodeGTF: gencodeGTF - gencodeTableBrowser: gencodeTableBrowser - repMaskBedFile: repMaskBEDFile - fileList1: fileListFile1 - - scatter: [repFamilySam, rmRepSam] - scatterMethod: dotproduct - - out: - - deduplicatedRmDupSam - - deduplicatedPreRmDupSam - - parsedFile - - doneFile - - step_concatenate_rmDup: - doc: "concatenates all rmdup sam files using cat" - run: concatenate.cwl - in: - files: step_deduplicate/deduplicatedRmDupSam - concatenated_output: - source: dataset - valueFrom: | - ${ - return self + ".rmDup.sam"; - } - out: - - concatenated - - step_concatenate_preRmDup: - doc: "concatenates all pre-rmduped sam files using cat" - run: concatenate.cwl - in: - files: step_deduplicate/deduplicatedPreRmDupSam - concatenated_output: - source: dataset - valueFrom: | - ${ - return self + ".preRmDup.sam"; - } - out: - - concatenated - - step_gzip_rmDup: - run: gzip.cwl - in: - input: step_concatenate_rmDup/concatenated - out: - - gzipped - - step_gzip_preRmDup: - run: gzip.cwl - in: - input: step_concatenate_preRmDup/concatenated - out: - - gzipped - - step_combine_parsed: - doc: "concatenates all final statistics using custom perl script" - run: combine.cwl - in: - files: step_deduplicate/parsedFile - outputFile: - source: dataset - valueFrom: | - ${ - return self + ".parsed"; - } - out: - - output diff --git a/cwl/wf_ecliprepmap_se.cwl b/cwl/wf_ecliprepmap_se.cwl deleted file mode 100755 index 1506207..0000000 --- a/cwl/wf_ecliprepmap_se.cwl +++ /dev/null @@ -1,170 +0,0 @@ -#!/usr/bin/env cwltool - -cwlVersion: v1.0 -class: Workflow - -requirements: - - class: SubworkflowFeatureRequirement - - class: MultipleInputFeatureRequirement - - class: InlineJavascriptRequirement - - class: StepInputExpressionRequirement - -inputs: - - dataset: - type: string - - barcode1r1FastqGz: - type: File - barcode1rmRepBam: - type: File - - barcode1Inputr1FastqGz: - type: File - barcode1InputrmRepBam: - type: File - - bowtie2_db: - type: Directory - bowtie2_prefix: - type: string - fileListFile1: - type: File - - gencodeGTF: - type: File - gencodeTableBrowser: - type: File - repMaskBEDFile: - type: File - - prefixes: - type: string[] - default: [ - "AA","AC","AG","AT","AN", - "CA","CC","CG","CT","CN", - "GA","GC","GG","GT","GN", - "TA","TC","TG","TT","TN", - "NA","NC","NG","NT","NN" - ] - se_or_pe: - type: string - default: SE - -outputs: - - - ### PRE RMDUPED SAM FILE FINAL OUTPUTS ### - - - output_ip_concatenated_pre_rmDup_sam_file: - type: File - outputSource: step_ecliprepmap_barcode1/output_concatenated_preRmDup_sam_file - output_input_concatenated_pre_rmDup_sam_file: - type: File - outputSource: step_ecliprepmap_input/output_concatenated_preRmDup_sam_file - - - ### RMDUPED SAM FILE FINAL OUTPUTS ### - - output_barcode1_concatenated_rmDup_sam_file: - type: File - outputSource: step_ecliprepmap_barcode1/output_concatenated_rmDup_sam_file - output_input_concatenated_rmDup_sam_file: - type: File - outputSource: step_ecliprepmap_input/output_concatenated_rmDup_sam_file - - - ### FINAL PARSED STATISTICS FILES ### - - output_ip_parsed: - type: File - outputSource: step_ecliprepmap_barcode1/output_combined_parsed_file - output_input_parsed: - type: File - outputSource: step_ecliprepmap_input/output_combined_parsed_file - output_nopipes: - type: File - outputSource: step_calculate_fold_change_from_parsed_files/out_file_nopipes_file - output_withpipes: - type: File - outputSource: step_calculate_fold_change_from_parsed_files/out_file_withpipes_file - -steps: - -########################################################################### -# Repeat-map single-end IP (1 sample) -########################################################################### - - step_ecliprepmap_barcode1: - run: wf_ecliprepmap_se_1barcode.cwl - in: - dataset: - source: barcode1r1FastqGz - valueFrom: | - ${ - return self.nameroot + ".barcode1"; - } - r1FastqGz: barcode1r1FastqGz - rmRepBam: barcode1rmRepBam - bowtie2_db: bowtie2_db - bowtie2_prefix: bowtie2_prefix - fileListFile1: fileListFile1 - gencodeGTF: gencodeGTF - gencodeTableBrowser: gencodeTableBrowser - repMaskBEDFile: repMaskBEDFile - prefixes: prefixes - se_or_pe: se_or_pe - out: - - output_repeat_mapped_sam_file - - output_rmDup_sam_files - - output_pre_rmDup_sam_files - - output_concatenated_rmDup_sam_file - - output_concatenated_preRmDup_sam_file - - output_parsed_files - - output_combined_parsed_file - -########################################################################### -# Repeat-map input sample (1 sample) -########################################################################### - - step_ecliprepmap_input: - run: wf_ecliprepmap_se_1barcode.cwl - in: - dataset: - source: barcode1Inputr1FastqGz - valueFrom: | - ${ - return self.nameroot + ".input"; - } - r1FastqGz: barcode1Inputr1FastqGz - bowtie2_db: bowtie2_db - bowtie2_prefix: bowtie2_prefix - fileListFile1: fileListFile1 - rmRepBam: barcode1InputrmRepBam - gencodeGTF: gencodeGTF - gencodeTableBrowser: gencodeTableBrowser - repMaskBEDFile: repMaskBEDFile - prefixes: prefixes - se_or_pe: se_or_pe - out: - - output_repeat_mapped_sam_file - - output_rmDup_sam_files - - output_pre_rmDup_sam_files - - output_concatenated_rmDup_sam_file - - output_concatenated_preRmDup_sam_file - - output_parsed_files - - output_combined_parsed_file - -########################################################################### -# Combine parsed files and calculate fold change/entropy -########################################################################### - - step_calculate_fold_change_from_parsed_files: - run: calculate_fold_change_from_parsed_files.cwl - in: - ip_parsed_file: step_ecliprepmap_barcode1/output_combined_parsed_file - input_parsed_file: step_ecliprepmap_input/output_combined_parsed_file - out: - - out_file_nopipes_file - - out_file_withpipes_file diff --git a/cwl/wf_ecliprepmap_se_1barcode.cwl b/cwl/wf_ecliprepmap_se_1barcode.cwl deleted file mode 100755 index ad37d49..0000000 --- a/cwl/wf_ecliprepmap_se_1barcode.cwl +++ /dev/null @@ -1,209 +0,0 @@ -#!/usr/bin/env cwltool - -cwlVersion: v1.0 -class: Workflow - -requirements: - - class: ScatterFeatureRequirement - - class: InlineJavascriptRequirement - - class: StepInputExpressionRequirement - -inputs: - - dataset: - type: string - - r1FastqGz: - doc: "Read1 after trimming inline barcodes" - type: File - rmRepBam: - doc: "BAM file after removing repetitive elements (after 2nd STAR mapping)" - type: File - - bowtie2_db: - type: Directory - bowtie2_prefix: - type: string - - fileListFile1: - type: File - - gencodeGTF: - type: File - gencodeTableBrowser: - type: File - repMaskBEDFile: - type: File - - prefixes: - type: string[] - default: [ - "AA","AC","AG","AT","AN", - "CA","CC","CG","CT","CN", - "GA","GC","GG","GT","GN", - "TA","TC","TG","TT","TN", - "NA","NC","NG","NT","NN" - ] - - se_or_pe: - type: string - default: SE - -outputs: - - output_repeat_mapped_sam_file: - doc: "first mapped SAM-like file to bowtie2" - type: File - outputSource: step_map_repetitive_elements/rep_sam - - output_rmDup_sam_files: - doc: "duplicate removed SAM-like files, one for each barcode" - type: File[] - outputSource: step_deduplicate/deduplicatedRmDupSam - - output_pre_rmDup_sam_files: - doc: "pre-duplicate removed SAM-like files, one for each barcode" - type: File[] - outputSource: step_deduplicate/deduplicatedPreRmDupSam - - output_concatenated_rmDup_sam_file: - doc: "Final remove duplicate SAM-like file" - type: File - outputSource: step_gzip_rmDup/gzipped - - output_concatenated_preRmDup_sam_file: - doc: "Final duplicated SAM-like file" - type: File - outputSource: step_gzip_preRmDup/gzipped - - output_parsed_files: - doc: "Output file containing read stats for each UMI prefix" - type: File[] - outputSource: step_deduplicate/parsedFile - - output_combined_parsed_file: - doc: "Combined output file containing read stats for all UMI prefixes" - type: File - outputSource: step_combine_parsed/output - -steps: - - step_map_repetitive_elements: - run: map_repetitive_elements_se.cwl - in: - read1: r1FastqGz - bowtie2_db: bowtie2_db - bowtie2_prefix: bowtie2_prefix - file_list_file: fileListFile1 - out: - - rep_sam - - step_splitbam_repsam: - run: splitbam.cwl - in: - sam_file: step_map_repetitive_elements/rep_sam - se_or_pe: se_or_pe - out: - - repsam_s - - step_splitbam_rmrepbam: - run: splitbam.cwl - in: - sam_file: rmRepBam - se_or_pe: se_or_pe - out: - - repsam_s - - step_getpair: - doc: | - Given a prefix (AA, AC, ... NN), return the rep and rmrep pairs - belonging to each prefix. Each file contains reads whose first 2nt - of its umi matches the prefix. - run: getpair.cwl - in: - rep_s: step_splitbam_repsam/repsam_s - rmrep_s: step_splitbam_rmrepbam/repsam_s - prefix: prefixes - scatter: prefix - out: - - prefixrep - - prefixrmrep - - step_deduplicate: - doc: | - Takes the repeat-mapped sam-like file, and the remove-replicate - bam file, and uses UMI information to remove PCR duplicates. - run: deduplicate.cwl - in: - repFamilySam: step_getpair/prefixrep - rmRepSam: step_getpair/prefixrmrep - se_or_pe: se_or_pe - gencodeGTF: gencodeGTF - gencodeTableBrowser: gencodeTableBrowser - repMaskBedFile: repMaskBEDFile - fileList1: fileListFile1 - - scatter: [repFamilySam, rmRepSam] - scatterMethod: dotproduct - - out: - - deduplicatedRmDupSam - - deduplicatedPreRmDupSam - - parsedFile - - doneFile - - step_concatenate_rmDup: - doc: "concatenates all rmdup sam files using cat" - run: concatenate.cwl - in: - files: step_deduplicate/deduplicatedRmDupSam - concatenated_output: - source: dataset - valueFrom: | - ${ - return self + ".rmDup.sam"; - } - out: - - concatenated - - step_concatenate_preRmDup: - doc: "concatenates all pre-rmduped sam files using cat" - run: concatenate.cwl - in: - files: step_deduplicate/deduplicatedPreRmDupSam - concatenated_output: - source: dataset - valueFrom: | - ${ - return self + ".preRmDup.sam"; - } - out: - - concatenated - - step_gzip_rmDup: - run: gzip.cwl - in: - input: step_concatenate_rmDup/concatenated - out: - - gzipped - - step_gzip_preRmDup: - run: gzip.cwl - in: - input: step_concatenate_preRmDup/concatenated - out: - - gzipped - - step_combine_parsed: - doc: "concatenates all final statistics using custom perl script" - run: combine.cwl - in: - files: step_deduplicate/parsedFile - outputFile: - source: dataset - valueFrom: | - ${ - return self + ".parsed"; - } - out: - - output diff --git a/docs/config_mapping.md b/docs/config_mapping.md new file mode 100644 index 0000000..09ec940 --- /dev/null +++ b/docs/config_mapping.md @@ -0,0 +1,49 @@ +# Config Mapping: CWL YAML to Snakemake config.yaml + +## Why they looked different +The previous `config/config.yaml` only covered mini validation fixtures for staged conversion tests. +It did not yet expose the full CWL-style runtime inputs. + +## Current mapping approach +`config/config.yaml` includes CWL-compatible key names at top level. +You can also pass an existing CWL-style job YAML via: + +```bash +snakemake --config cwl_input_yaml=/path/to/job.yaml ... +``` + +When `cwl_input_yaml` is provided, mapped keys from that file override defaults in `config/config.yaml`. + +CWL YAML key -> Snakemake key: +- `dataset` -> `dataset` +- `barcode1r1FastqGz` -> `barcode1r1FastqGz` +- `barcode1rmRepBam` -> `barcode1rmRepBam` +- `barcode1Inputr1FastqGz` -> `barcode1Inputr1FastqGz` +- `barcode1InputrmRepBam` -> `barcode1InputrmRepBam` +- `bowtie2_db` -> `bowtie2_db` +- `bowtie2_prefix` -> `bowtie2_prefix` +- `fileListFile1` -> `fileListFile1` +- `gencodeGTF` -> `gencodeGTF` +- `gencodeTableBrowser` -> `gencodeTableBrowser` +- `repMaskBEDFile` -> `repMaskBEDFile` +- `prefixes` -> `prefixes` +- `se_or_pe` -> `se_or_pe` + +## Validation-only section +The `mini_validation` subsection is separate and only used for deterministic fixture testing of currently implemented stages. + +## `class/path` object handling +For any supported key with a value like: + +```yaml +someKey: + class: File + path: relative/or/absolute/path +``` + +the adapter uses `path` as the runtime value. +Relative paths are resolved relative to the job YAML file directory. + +## Practical note +For production runs, populate the top-level CWL-compatible keys with real file paths. +For conversion regression tests, keep using `mini_validation`. diff --git a/docs/data_inventory.md b/docs/data_inventory.md new file mode 100644 index 0000000..9a7dd9b --- /dev/null +++ b/docs/data_inventory.md @@ -0,0 +1,62 @@ +# Data Inventory (PPA1 Example + Required References) + +## Location and path semantics +- In this workspace, `data/` is a symlink to `/Volumes/X9Pro/Yeo/repetitive-element-pipeline`. +- All paths below are documented using repository-relative paths for reproducibility. + +## PPA1 run metadata +- Input YAML: + - `data/PPA1_rep1.yaml` + - `data/PPA1_rep1/REPELEMENTMAPPING_PPA1_rep1_INPUT.yaml` +- Workflow run metadata: + - `data/PPA1_rep1/REPELEMENTMAPPING_PPA1_rep1_OUTPUT.json` + - `data/PPA1_rep1/REPELEMENTMAPPING_PPA1_rep1_LOG.txt` + - `data/PPA1_rep1/REPELEMENTMAPPING_PPA1_rep1_WORKFLOW-wf_ecliprepmap_se` + - `data/PPA1_rep1/REPELEMENTMAPPING_PPA1_rep1_VERSION-0.1.0` + +## Required SE reference inputs + +| Path | Purpose | Size (bytes) | +|---|---|---:| +| `data/bowtie_reference/MASTER_filelist.wrepbaseandtRNA.fa.fixed.fa.UpdatedSimpleRepeat.{1,2,3,4,rev.1,rev.2}.bt2` | Bowtie2 index shards | varied | +| `data/MASTER_filelist.wrepbaseandtRNA.enst2id.fixed.UpdatedSimpleRepeat` | ENST-to-family/type mapping list | 616168 | +| `data/gencode.v19.chr_patch_hapl_scaff.annotation.gtf` | Gencode annotation | 1199970581 | +| `data/gencode.v19.chr_patch_hapl_scaff.annotation.gtf.parsed_ucsc_tableformat` | Parsed Gencode table format | 41595212 | +| `data/RepeatMask.bed` | RepeatMasker regions | 180955350 | + +Additional provided references (not required for minimal SE CWL execution path): +- `data/ALLRepBase_elements.id_table.FULL` +- `data/genelists.chrM.wchr.txt` +- `data/mirbase.v20.hg19.gff3` + +## PPA1 example result artifacts +`data/PPA1_rep1/results/` contains expected compressed outputs from the SE workflow, including: +- Parsed summaries: `*.parsed.gz`, `*.reparsed.tsv.gz` +- Dedup outputs: `*.preRmDup.sam.gz`, `*.rmDup.sam.gz` +- Fold-change tables: `*.nopipes.tsv.gz`, `*.withpipes.tsv.gz` +- mpileup-short files: `*.mpileup.short.gz` + +A full expected file fingerprint table is recorded in: +- `docs/ppa1_expected_outputs_manifest.tsv` + +## Observed line-level structures (sampled) + +### SAM-like rmDup output (`*.rmDup.sam.gz`) +- Contains standard SAM columns plus appended annotations such as `RepFamily`/`UniqueGenomic` and normalized family labels. + +### Parsed summary (`*.parsed.gz`) +- Starts with `#READINFO` header lines, e.g. `AllReads`, `UsableReads`, `GenomicReads`, `RepFamilyReads`. + +### Reparsed fold-change input (`*.reparsed.nopipes.tsv.gz`) +- Header columns: + - `element`, `IP_read_num`, `IP_clip_rpr`, `Input_read_num`, `Input_clip_rpr`, `Fold_enrichment`, `Information_content` + +## Reproducible manifest generation +Command: +```bash +python3 scripts/generate_ppa1_manifest.py \ + --input-dir data/PPA1_rep1/results \ + --output docs/ppa1_expected_outputs_manifest.tsv +``` + +This records `relative_path`, `size_bytes`, and `sha256` for each expected output file. diff --git a/docs/mini_fixtures.md b/docs/mini_fixtures.md new file mode 100644 index 0000000..8244cd6 --- /dev/null +++ b/docs/mini_fixtures.md @@ -0,0 +1,33 @@ +# Mini Fixture Set + +## Purpose +These mini fixtures provide deterministic, small test data derived from real `PPA1_rep1` outputs so we can validate Perl-to-Python and CWL-to-Snakemake conversion steps quickly. + +## Layout +- Source mini inputs: + - `tests/fixtures/mini/source/ip.preRmDup.sam.mini.gz` + - `tests/fixtures/mini/source/input.preRmDup.sam.mini.gz` +- Expected mini outputs: + - `tests/fixtures/mini/expected/ip.rmDup.sam.mini.gz` + - `tests/fixtures/mini/expected/input.rmDup.sam.mini.gz` + - `tests/fixtures/mini/expected/ip.parsed.mini.gz` + - `tests/fixtures/mini/expected/input.parsed.mini.gz` + - `tests/fixtures/mini/expected/ip.reparsed.nopipes.mini.tsv.gz` + - `tests/fixtures/mini/expected/ip.reparsed.withpipes.mini.tsv.gz` +- Checksum manifest: + - `tests/fixtures/mini/manifest.tsv` + +## Generation +```bash +python3 scripts/make_mini_fixtures.py +``` + +Defaults: +- SAM-like files: first 3000 lines +- Parsed files: all `#READINFO` + first 60 `TOTAL` + first 200 `ELEMENT` +- TSV files: header + first 300 data rows + +## Validation +```bash +python3 -m pytest -q tests/test_mini_fixtures_manifest.py +``` diff --git a/docs/output_format_spec.md b/docs/output_format_spec.md new file mode 100644 index 0000000..25415be --- /dev/null +++ b/docs/output_format_spec.md @@ -0,0 +1,75 @@ +# Output Format Specification (Observed + Code-Aligned) + +## Source basis +- Observed files under `data/PPA1_rep1/results/*.gz`. +- CWL tool wiring and Perl/Python scripts. + +## 1) `.parsed` and `.reparsed.tsv`-style summary files + +### Header lines +Common readinfo-like rows begin with `#READINFO`. +Observed keys: +- `AllReads` (in `.parsed`, not always present in reparsed files) +- `UsableReads` +- `GenomicReads` +- `RepFamilyReads` + +### Body rows +Two row classes: +- `TOTAL\t\t\t` +- `ELEMENT\t\t\t\t\t` + +Notes: +- Some files labeled `.reparsed.tsv.gz` retain `TOTAL` rows but may omit initial `AllReads` line. +- Floating precision differs by script/formatting (`sprintf` in Perl; pandas/float formatting in Python). + +## 2) `.nopipes.tsv` and `.withpipes.tsv` +Tab-separated with header: +- `element` +- `IP_read_num` +- `IP_clip_rpr` +- `Input_read_num` +- `Input_clip_rpr` +- `Fold_enrichment` +- `Information_content` + +Meaning: +- `clip_rpr` terms are read proportion metrics. +- `Fold_enrichment = IP_clip_rpr / Input_clip_rpr` (with pseudocount handling in upstream script). +- `Information_content = IP_clip_rpr * log2(IP_clip_rpr / Input_clip_rpr)`. +- `.nopipes.tsv` excludes ambiguous multi-family mappings (contains no `|` in element family labels). +- `.withpipes.tsv` includes both ambiguous and unambiguous families. + +## 3) SAM-like pre-rmdup and rmdup outputs +Examples: +- `*.preRmDup.sam.gz` +- `*.rmDup.sam.gz` + +Characteristics: +- Tab-delimited SAM fields preserved. +- Additional right-side annotations are appended beyond canonical SAM tags. +- `ZZ:Z:` tags may include pipe-delimited ENST sets. +- rmdup files append classification fields such as: + - `RepFamily` / `UniqueGenomic` + - normalized family/element label string + +## 4) mpileup-short outputs +Examples: +- `*.rmDup.sam.gz.tmp.RNA18S.bam.sorted.bam.mpileup.short.gz` +- `*.rmDup.sam.gz.tmp.RNA28S.bam.sorted.bam.mpileup.short.gz` + +Observed columns: +1. transcript/accession id +2. 1-based position +3. reference base +4. depth/count + +## 5) Determinism/tolerance for validation +Validation should be tolerant for: +- floating point formatting differences +- row ordering differences where hash/dict iteration affects output order + +Validation should be strict for: +- schema/column set +- primary key identity (element/family labels) +- integer read counts diff --git a/docs/ppa1_expected_outputs_manifest.tsv b/docs/ppa1_expected_outputs_manifest.tsv new file mode 100644 index 0000000..f637b6f --- /dev/null +++ b/docs/ppa1_expected_outputs_manifest.tsv @@ -0,0 +1,19 @@ +relative_path size_bytes sha256 +data/PPA1_rep1/results/MP2.PPA1_IN1.umi.r1.fqTrTr.sorted.fq.input.parsed.gz 65365973 6b7dd5f40ad086e407379e757a508e9925d23b7de47c12198b9cde2377618678 +data/PPA1_rep1/results/MP2.PPA1_IN1.umi.r1.fqTrTr.sorted.fq.input.preRmDup.sam.gz 2962472583 ef0ffdfdf080afd5b1b6a1c293229733b2e141fbd6ae99965d94da0bbde65867 +data/PPA1_rep1/results/MP2.PPA1_IN1.umi.r1.fqTrTr.sorted.fq.input.rmDup.sam.gz 1799892699 bc3ab8725324af85fe1e59cdf97bb464b3b24ce5814886823de0ceab0d087e98 +data/PPA1_rep1/results/MP2.PPA1_IN1.umi.r1.fqTrTr.sorted.fq.input.rmDup.sam.gz.tmp.RNA18S.bam.sorted.bam.mpileup.short.gz 12786 d424ef50cd119611a6d5c2ebcda290aef239d008aa512ab84516a71f8d68a9fd +data/PPA1_rep1/results/MP2.PPA1_IN1.umi.r1.fqTrTr.sorted.fq.input.rmDup.sam.gz.tmp.RNA28S.bam.sorted.bam.mpileup.short.gz 34240 177a27117d12746c29d0962f0d7c2ced84a94956c5a66e609e886260b93afd08 +data/PPA1_rep1/results/MP2.PPA1_IN1.umi.r1.fqTrTr.sorted.fq.input.rmDup.sam.gz.tmp.rRNA_extra.bam.sorted.bam.mpileup.short.gz 38440 08dae7e0e1bf4a80f26db4c509c94d5f3e694860b9e75da08b6a9e03a11c57ea +data/PPA1_rep1/results/MP2.PPA1_IN1.umi.r1.fqTrTr.sorted.fq.input.rmDup.sam.reparsed.tsv.gz 2984352 5b272f3bac56b2817e0d14b69e80a8576ad6c401ff9dcbcedc18db5554c071fb +data/PPA1_rep1/results/MP2.PPA1_IP1.umi.r1.fqTrTr.sorted.fq.barcode1.nopipes.tsv.gz 8008 49059581de39e8aa90f1dfcfee37fce77f6bead55dd6bba6dcdf33db1b09dae1 +data/PPA1_rep1/results/MP2.PPA1_IP1.umi.r1.fqTrTr.sorted.fq.barcode1.parsed.gz 29182643 ad30c5775ef021e33fa6a34f3ab22fa6ebcee0d5bc8a272182b2aa3d644c6c65 +data/PPA1_rep1/results/MP2.PPA1_IP1.umi.r1.fqTrTr.sorted.fq.barcode1.preRmDup.sam.gz 2192050419 48c0b76ae20cd3a16797eb3f47f81cd0b9612df356fe2e3fd412ea5c98536ee2 +data/PPA1_rep1/results/MP2.PPA1_IP1.umi.r1.fqTrTr.sorted.fq.barcode1.rmDup.sam.gz 1288407573 59a823800ba37d469494db7c8b92ca54ecdc7dab14bdaa64b5f7548ff42fe041 +data/PPA1_rep1/results/MP2.PPA1_IP1.umi.r1.fqTrTr.sorted.fq.barcode1.rmDup.sam.gz.tmp.RNA18S.bam.sorted.bam.mpileup.short.gz 12587 69c387d959f8df7e1b2a91a77d8d916985b6b43cb34053a06c7dc93812662383 +data/PPA1_rep1/results/MP2.PPA1_IP1.umi.r1.fqTrTr.sorted.fq.barcode1.rmDup.sam.gz.tmp.RNA28S.bam.sorted.bam.mpileup.short.gz 32816 ea8e60c87cc811af71154c8d5ac8dcf23abf780c2e2f66f3f670af6659629354 +data/PPA1_rep1/results/MP2.PPA1_IP1.umi.r1.fqTrTr.sorted.fq.barcode1.rmDup.sam.gz.tmp.rRNA_extra.bam.sorted.bam.mpileup.short.gz 37907 406de4534758efc703af97939d2b702d3f9a44f09d69a8f81930a5b00109eef5 +data/PPA1_rep1/results/MP2.PPA1_IP1.umi.r1.fqTrTr.sorted.fq.barcode1.rmDup.sam.reparsed.nopipes.tsv.gz 7236 b67e5b412cf433fb94be392ca45428253d549fae82f4b481e248c7d0b396dbe5 +data/PPA1_rep1/results/MP2.PPA1_IP1.umi.r1.fqTrTr.sorted.fq.barcode1.rmDup.sam.reparsed.tsv.gz 3661993 d102231a4c531ac0a8e952d135963e71f953e9afcab5acef4ba24d129a13db7b +data/PPA1_rep1/results/MP2.PPA1_IP1.umi.r1.fqTrTr.sorted.fq.barcode1.rmDup.sam.reparsed.withpipes.tsv.gz 42520 eaf938eb4491f19bf898e52497431def45aaee4fed84d8fcf8bd363f995d7852 +data/PPA1_rep1/results/MP2.PPA1_IP1.umi.r1.fqTrTr.sorted.fq.barcode1.withpipes.tsv.gz 44057 3cc5817f355a8addc1627500f61fd2b220f6790ff89f7905f8ed25a66d3dfe37 diff --git a/docs/snakemake_se_foundation.md b/docs/snakemake_se_foundation.md new file mode 100644 index 0000000..dce065e --- /dev/null +++ b/docs/snakemake_se_foundation.md @@ -0,0 +1,28 @@ +# Snakemake SE Foundation + +## Implemented stages +1. Split SAM by UMI prefix (`split_bam_to_subfiles_SEorPE.py`). +2. Merge parsed statistics files (`merge_multiple_parsed_files.simplified_20191022.py`). + +## Files +- `Snakefile` +- `config/config.yaml` +- `workflow/rules/se_foundation.smk` +- `workflow/scripts/verify_split_manifest.py` +- `workflow/scripts/verify_merge_against_expected.py` + +## Validation strategy +- Split stage output files are validated against + `tests/fixtures/mini/expected/split.expected.manifest.tsv`. +- Merge stage output is validated against + `tests/fixtures/mini/expected/merged.expected.parsed`. + +## Run +```bash +HOME=$(pwd) ./.conda-env/bin/snakemake -j1 -p -F all +``` + +## Next extensions +- Add mapping stage rules. +- Add deduplication stage rules. +- Add full end-to-end stage chaining for production data. diff --git a/documentation/lab_meeting_slides_rep_element_pipeline-20171018.pdf b/documentation/lab_meeting_slides_rep_element_pipeline-20171018.pdf deleted file mode 100644 index 1259741..0000000 Binary files a/documentation/lab_meeting_slides_rep_element_pipeline-20171018.pdf and /dev/null differ diff --git a/documentation/lab_meeting_slides_rep_element_pipeline.pdf b/documentation/lab_meeting_slides_rep_element_pipeline.pdf deleted file mode 100644 index 1259741..0000000 Binary files a/documentation/lab_meeting_slides_rep_element_pipeline.pdf and /dev/null differ diff --git a/scripts/generate_ppa1_manifest.py b/scripts/generate_ppa1_manifest.py new file mode 100755 index 0000000..d8a4841 --- /dev/null +++ b/scripts/generate_ppa1_manifest.py @@ -0,0 +1,39 @@ +#!/usr/bin/env python3 +from __future__ import annotations + +import argparse +import hashlib +from pathlib import Path + + +def sha256_file(path: Path, chunk_size: int = 1024 * 1024) -> str: + h = hashlib.sha256() + with path.open('rb') as fh: + while True: + data = fh.read(chunk_size) + if not data: + break + h.update(data) + return h.hexdigest() + + +def main() -> None: + parser = argparse.ArgumentParser(description='Generate expected output manifest for PPA1 example files') + parser.add_argument('--input-dir', default='data/PPA1_rep1/results', help='Directory containing expected result files') + parser.add_argument('--output', default='docs/ppa1_expected_outputs_manifest.tsv', help='TSV output path') + args = parser.parse_args() + + in_dir = Path(args.input_dir) + out = Path(args.output) + + files = sorted(p for p in in_dir.iterdir() if p.is_file()) + out.parent.mkdir(parents=True, exist_ok=True) + + with out.open('w', encoding='utf-8') as fh: + fh.write('relative_path\tsize_bytes\tsha256\n') + for path in files: + fh.write(f"{path.as_posix()}\t{path.stat().st_size}\t{sha256_file(path)}\n") + + +if __name__ == '__main__': + main() diff --git a/scripts/generate_python_expected_artifacts.py b/scripts/generate_python_expected_artifacts.py new file mode 100755 index 0000000..16ac2c2 --- /dev/null +++ b/scripts/generate_python_expected_artifacts.py @@ -0,0 +1,72 @@ +#!/usr/bin/env python3 +from __future__ import annotations + +import csv +import gzip +import hashlib +import shutil +import subprocess +from pathlib import Path + + +REPO = Path(__file__).resolve().parents[1] +SOURCE_GZ = REPO / 'tests/fixtures/mini/source/ip.preRmDup.sam.mini.gz' +MERGE_IN1 = REPO / 'tests/fixtures/mini/source/merge_input_1.parsed_v2.txt' +MERGE_IN2 = REPO / 'tests/fixtures/mini/source/merge_input_2.parsed_v2.txt' +SPLIT_MANIFEST = REPO / 'tests/fixtures/mini/expected/split.expected.manifest.tsv' +MERGED_EXPECTED = REPO / 'tests/fixtures/mini/expected/merged.expected.parsed' + + +def sha256_file(path: Path) -> str: + h = hashlib.sha256() + with path.open('rb') as fh: + while True: + chunk = fh.read(1024 * 1024) + if not chunk: + break + h.update(chunk) + return h.hexdigest() + + +def main() -> None: + tmp = REPO / 'results/.expected_tmp' + if tmp.exists(): + shutil.rmtree(tmp) + tmp.mkdir(parents=True) + + sam = tmp / 'ip.preRmDup.sam.mini.sam' + with gzip.open(SOURCE_GZ, 'rt', encoding='utf-8', errors='replace') as src, sam.open('w', encoding='utf-8') as out: + out.write(src.read()) + + split_dir = tmp / 'split' + split_dir.mkdir() + subprocess.run( + ['python3', str(REPO / 'bin/python/split_bam_to_subfiles_SEorPE.py'), str(sam), 'SE'], + cwd=split_dir, + check=True, + ) + + SPLIT_MANIFEST.parent.mkdir(parents=True, exist_ok=True) + with SPLIT_MANIFEST.open('w', encoding='utf-8', newline='') as fh: + writer = csv.writer(fh, delimiter='\t') + writer.writerow(['filename', 'size_bytes', 'sha256']) + for path in sorted(split_dir.glob('*.tmp')): + writer.writerow([path.name, path.stat().st_size, sha256_file(path)]) + + MERGED_EXPECTED.parent.mkdir(parents=True, exist_ok=True) + subprocess.run( + [ + 'python3', + str(REPO / 'bin/python/merge_multiple_parsed_files.simplified_20191022.py'), + str(MERGED_EXPECTED), + str(MERGE_IN1), + str(MERGE_IN2), + ], + check=True, + ) + + shutil.rmtree(tmp) + + +if __name__ == '__main__': + main() diff --git a/scripts/make_mini_fixtures.py b/scripts/make_mini_fixtures.py new file mode 100755 index 0000000..3d03809 --- /dev/null +++ b/scripts/make_mini_fixtures.py @@ -0,0 +1,126 @@ +#!/usr/bin/env python3 +from __future__ import annotations + +import argparse +import gzip +import hashlib +from pathlib import Path +from typing import Iterable + + +def iter_lines_gz(path: Path) -> Iterable[str]: + with gzip.open(path, 'rt', encoding='utf-8', errors='replace') as fh: + for line in fh: + yield line + + +def write_gz(path: Path, lines: Iterable[str]) -> None: + path.parent.mkdir(parents=True, exist_ok=True) + with gzip.open(path, 'wt', encoding='utf-8') as out: + for line in lines: + out.write(line) + + +def head_lines(lines: Iterable[str], n: int) -> list[str]: + out: list[str] = [] + for line in lines: + out.append(line) + if len(out) >= n: + break + return out + + +def parsed_subset(lines: Iterable[str], total_rows: int, element_rows: int) -> list[str]: + out: list[str] = [] + total_seen = 0 + element_seen = 0 + for line in lines: + if line.startswith('#READINFO'): + out.append(line) + continue + if line.startswith('TOTAL\t') and total_seen < total_rows: + out.append(line) + total_seen += 1 + continue + if line.startswith('ELEMENT\t') and element_seen < element_rows: + out.append(line) + element_seen += 1 + continue + if total_seen >= total_rows and element_seen >= element_rows: + break + return out + + +def tsv_subset(lines: Iterable[str], n_data_rows: int) -> list[str]: + out: list[str] = [] + for line in lines: + out.append(line) + if len(out) == 1: + continue + if len(out) - 1 >= n_data_rows: + break + return out + + +def sha256_file(path: Path) -> str: + h = hashlib.sha256() + with path.open('rb') as fh: + while True: + data = fh.read(1024 * 1024) + if not data: + break + h.update(data) + return h.hexdigest() + + +def main() -> None: + p = argparse.ArgumentParser(description='Create mini deterministic fixtures from PPA1 example outputs') + p.add_argument('--source-dir', default='data/PPA1_rep1/results') + p.add_argument('--out-dir', default='tests/fixtures/mini') + p.add_argument('--sam-lines', type=int, default=3000) + p.add_argument('--parsed-total-rows', type=int, default=60) + p.add_argument('--parsed-element-rows', type=int, default=200) + p.add_argument('--tsv-rows', type=int, default=300) + args = p.parse_args() + + source = Path(args.source_dir) + out = Path(args.out_dir) + src_out = out / 'source' + expected_out = out / 'expected' + + mapping = { + 'MP2.PPA1_IP1.umi.r1.fqTrTr.sorted.fq.barcode1.preRmDup.sam.gz': ('source/ip.preRmDup.sam.mini.gz', 'sam'), + 'MP2.PPA1_IN1.umi.r1.fqTrTr.sorted.fq.input.preRmDup.sam.gz': ('source/input.preRmDup.sam.mini.gz', 'sam'), + 'MP2.PPA1_IP1.umi.r1.fqTrTr.sorted.fq.barcode1.rmDup.sam.gz': ('expected/ip.rmDup.sam.mini.gz', 'sam'), + 'MP2.PPA1_IN1.umi.r1.fqTrTr.sorted.fq.input.rmDup.sam.gz': ('expected/input.rmDup.sam.mini.gz', 'sam'), + 'MP2.PPA1_IP1.umi.r1.fqTrTr.sorted.fq.barcode1.parsed.gz': ('expected/ip.parsed.mini.gz', 'parsed'), + 'MP2.PPA1_IN1.umi.r1.fqTrTr.sorted.fq.input.parsed.gz': ('expected/input.parsed.mini.gz', 'parsed'), + 'MP2.PPA1_IP1.umi.r1.fqTrTr.sorted.fq.barcode1.rmDup.sam.reparsed.nopipes.tsv.gz': ('expected/ip.reparsed.nopipes.mini.tsv.gz', 'tsv'), + 'MP2.PPA1_IP1.umi.r1.fqTrTr.sorted.fq.barcode1.rmDup.sam.reparsed.withpipes.tsv.gz': ('expected/ip.reparsed.withpipes.mini.tsv.gz', 'tsv'), + } + + generated: list[Path] = [] + for src_name, (rel_dst, kind) in mapping.items(): + src = source / src_name + dst = out / rel_dst + lines = iter_lines_gz(src) + if kind == 'sam': + subset = head_lines(lines, args.sam_lines) + elif kind == 'parsed': + subset = parsed_subset(lines, args.parsed_total_rows, args.parsed_element_rows) + elif kind == 'tsv': + subset = tsv_subset(lines, args.tsv_rows) + else: + raise ValueError(kind) + write_gz(dst, subset) + generated.append(dst) + + manifest = out / 'manifest.tsv' + with manifest.open('w', encoding='utf-8') as fh: + fh.write('relative_path\tsize_bytes\tsha256\n') + for path in sorted(generated): + fh.write(f"{path.as_posix()}\t{path.stat().st_size}\t{sha256_file(path)}\n") + + +if __name__ == '__main__': + main() diff --git a/tests/fixtures/mini/cwl_job_example.yaml b/tests/fixtures/mini/cwl_job_example.yaml new file mode 100644 index 0000000..abfb70b --- /dev/null +++ b/tests/fixtures/mini/cwl_job_example.yaml @@ -0,0 +1,38 @@ +dataset: fixture_dataset + +barcode1r1FastqGz: + class: File + path: source/ip.preRmDup.sam.mini.gz + +barcode1rmRepBam: + class: File + path: /abs/path/ip.bam + +barcode1Inputr1FastqGz: + class: File + path: source/input.preRmDup.sam.mini.gz + +barcode1InputrmRepBam: + class: File + path: /abs/path/input.bam + +bowtie2_db: + class: Directory + path: refs/bowtie2_index + +bowtie2_prefix: MASTER_FILELIST +fileListFile1: + class: File + path: refs/filelist.tsv +gencodeGTF: + class: File + path: refs/gencode.gtf +gencodeTableBrowser: + class: File + path: refs/gencode.table.tsv +repMaskBEDFile: + class: File + path: refs/repmask.bed + +prefixes: [AA, AC, TT] +se_or_pe: SE diff --git a/tests/fixtures/mini/expected/input.parsed.mini.gz b/tests/fixtures/mini/expected/input.parsed.mini.gz new file mode 100644 index 0000000..c1de19d Binary files /dev/null and b/tests/fixtures/mini/expected/input.parsed.mini.gz differ diff --git a/tests/fixtures/mini/expected/input.rmDup.sam.mini.gz b/tests/fixtures/mini/expected/input.rmDup.sam.mini.gz new file mode 100644 index 0000000..9e41776 Binary files /dev/null and b/tests/fixtures/mini/expected/input.rmDup.sam.mini.gz differ diff --git a/tests/fixtures/mini/expected/ip.parsed.mini.gz b/tests/fixtures/mini/expected/ip.parsed.mini.gz new file mode 100644 index 0000000..f0794d4 Binary files /dev/null and b/tests/fixtures/mini/expected/ip.parsed.mini.gz differ diff --git a/tests/fixtures/mini/expected/ip.reparsed.nopipes.mini.tsv.gz b/tests/fixtures/mini/expected/ip.reparsed.nopipes.mini.tsv.gz new file mode 100644 index 0000000..ffdd9bf Binary files /dev/null and b/tests/fixtures/mini/expected/ip.reparsed.nopipes.mini.tsv.gz differ diff --git a/tests/fixtures/mini/expected/ip.reparsed.withpipes.mini.tsv.gz b/tests/fixtures/mini/expected/ip.reparsed.withpipes.mini.tsv.gz new file mode 100644 index 0000000..609a915 Binary files /dev/null and b/tests/fixtures/mini/expected/ip.reparsed.withpipes.mini.tsv.gz differ diff --git a/tests/fixtures/mini/expected/ip.rmDup.sam.mini.gz b/tests/fixtures/mini/expected/ip.rmDup.sam.mini.gz new file mode 100644 index 0000000..a26e7a0 Binary files /dev/null and b/tests/fixtures/mini/expected/ip.rmDup.sam.mini.gz differ diff --git a/tests/fixtures/mini/expected/merged.expected.parsed b/tests/fixtures/mini/expected/merged.expected.parsed new file mode 100644 index 0000000..6f389ba --- /dev/null +++ b/tests/fixtures/mini/expected/merged.expected.parsed @@ -0,0 +1,8 @@ +#READINFO AllReads 300 +#READINFO UsableReads 240 0.8 +#READINFO GenomicReads 90 0.375 +#READINFO RepFamilyReads 150 0.625 +TOTAL RNA28S 90 0.375 +TOTAL RNA18S 60 0.25 +ELEMENT RNA28S 53 0.22083333333333333 RNA28S||ENST2 GENE2 +ELEMENT RNA18S 37 0.15416666666666667 RNA18S||ENST1 GENE1 diff --git a/tests/fixtures/mini/expected/split.expected.manifest.tsv b/tests/fixtures/mini/expected/split.expected.manifest.tsv new file mode 100644 index 0000000..0d47051 --- /dev/null +++ b/tests/fixtures/mini/expected/split.expected.manifest.tsv @@ -0,0 +1,26 @@ +filename size_bytes sha256 +AA.ip.preRmDup.sam.mini.sam.tmp 1287381 c30ea41eccd74859e6283d8fb54905829d3a0d0758bea939cee7357a00c7064f +AC.ip.preRmDup.sam.mini.sam.tmp 0 e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855 +AG.ip.preRmDup.sam.mini.sam.tmp 0 e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855 +AN.ip.preRmDup.sam.mini.sam.tmp 0 e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855 +AT.ip.preRmDup.sam.mini.sam.tmp 0 e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855 +CA.ip.preRmDup.sam.mini.sam.tmp 0 e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855 +CC.ip.preRmDup.sam.mini.sam.tmp 0 e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855 +CG.ip.preRmDup.sam.mini.sam.tmp 0 e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855 +CN.ip.preRmDup.sam.mini.sam.tmp 0 e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855 +CT.ip.preRmDup.sam.mini.sam.tmp 0 e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855 +GA.ip.preRmDup.sam.mini.sam.tmp 0 e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855 +GC.ip.preRmDup.sam.mini.sam.tmp 0 e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855 +GG.ip.preRmDup.sam.mini.sam.tmp 0 e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855 +GN.ip.preRmDup.sam.mini.sam.tmp 0 e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855 +GT.ip.preRmDup.sam.mini.sam.tmp 0 e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855 +NA.ip.preRmDup.sam.mini.sam.tmp 0 e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855 +NC.ip.preRmDup.sam.mini.sam.tmp 0 e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855 +NG.ip.preRmDup.sam.mini.sam.tmp 0 e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855 +NN.ip.preRmDup.sam.mini.sam.tmp 0 e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855 +NT.ip.preRmDup.sam.mini.sam.tmp 0 e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855 +TA.ip.preRmDup.sam.mini.sam.tmp 0 e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855 +TC.ip.preRmDup.sam.mini.sam.tmp 0 e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855 +TG.ip.preRmDup.sam.mini.sam.tmp 0 e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855 +TN.ip.preRmDup.sam.mini.sam.tmp 0 e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855 +TT.ip.preRmDup.sam.mini.sam.tmp 0 e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855 diff --git a/tests/fixtures/mini/manifest.tsv b/tests/fixtures/mini/manifest.tsv new file mode 100644 index 0000000..d424f47 --- /dev/null +++ b/tests/fixtures/mini/manifest.tsv @@ -0,0 +1,9 @@ +relative_path size_bytes sha256 +tests/fixtures/mini/expected/input.parsed.mini.gz 6492 63a47234eff4e54612065eb4386bb24c2a817165884c15b3aebb5245e199148d +tests/fixtures/mini/expected/input.rmDup.sam.mini.gz 174824 479a3774b8867e41292aff93a23c539cd220549522a2244a49e1b6d509f7f4a1 +tests/fixtures/mini/expected/ip.parsed.mini.gz 6292 fdc33e2b6f6a6b7df23f21c5bbffc87bfcd0e682d1309e541d703edf9e5c74dd +tests/fixtures/mini/expected/ip.reparsed.nopipes.mini.tsv.gz 7188 ab76d446267a23dad0f5d815ef0db2bc8d152155cbd7ac5d6ecf0be089ec5c72 +tests/fixtures/mini/expected/ip.reparsed.withpipes.mini.tsv.gz 12851 db976f71def48a10632aeaaa71226a6c19260c254da85c7f9e0afb0de3339255 +tests/fixtures/mini/expected/ip.rmDup.sam.mini.gz 185900 920ce5732417fcb46d9fbcf3c348e2dbf99f4d8fd74fe8c3248231f9855d77e8 +tests/fixtures/mini/source/input.preRmDup.sam.mini.gz 198136 24f3307ebfce353325bdea0813a54b2201ef83da12059af8aa3292a4fe705f2b +tests/fixtures/mini/source/ip.preRmDup.sam.mini.gz 204694 f6749227444bc327255cd98fe7e54be5c7c19daaacc90dd15ec2fbb996f9f9e6 diff --git a/tests/fixtures/mini/source/input.preRmDup.sam.mini.gz b/tests/fixtures/mini/source/input.preRmDup.sam.mini.gz new file mode 100644 index 0000000..eafe17a Binary files /dev/null and b/tests/fixtures/mini/source/input.preRmDup.sam.mini.gz differ diff --git a/tests/fixtures/mini/source/ip.preRmDup.sam.mini.gz b/tests/fixtures/mini/source/ip.preRmDup.sam.mini.gz new file mode 100644 index 0000000..f50c874 Binary files /dev/null and b/tests/fixtures/mini/source/ip.preRmDup.sam.mini.gz differ diff --git a/tests/fixtures/mini/source/merge_input_1.parsed_v2.txt b/tests/fixtures/mini/source/merge_input_1.parsed_v2.txt new file mode 100644 index 0000000..e1f91bf --- /dev/null +++ b/tests/fixtures/mini/source/merge_input_1.parsed_v2.txt @@ -0,0 +1,8 @@ +#READINFO All reads: 100 PCR duplicates removed: 20 Usable Remaining: 80 Usable from genomic mapping: 30 Usable from family mapping: 50 +#READINFO UsableReads 80 +#READINFO GenomicReads 30 0.375 +#READINFO RepFamilyReads 50 0.625 +TOTAL RNA18S 20 250000 +TOTAL RNA28S 30 375000 +ELEMENT RNA18S 12 150000 RNA18S||ENST1 GENE1 +ELEMENT RNA28S 18 225000 RNA28S||ENST2 GENE2 diff --git a/tests/fixtures/mini/source/merge_input_2.parsed_v2.txt b/tests/fixtures/mini/source/merge_input_2.parsed_v2.txt new file mode 100644 index 0000000..5863224 --- /dev/null +++ b/tests/fixtures/mini/source/merge_input_2.parsed_v2.txt @@ -0,0 +1,8 @@ +#READINFO All reads: 200 PCR duplicates removed: 40 Usable Remaining: 160 Usable from genomic mapping: 60 Usable from family mapping: 100 +#READINFO UsableReads 160 +#READINFO GenomicReads 60 0.375 +#READINFO RepFamilyReads 100 0.625 +TOTAL RNA18S 40 250000 +TOTAL RNA28S 60 375000 +ELEMENT RNA18S 25 156250 RNA18S||ENST1 GENE1 +ELEMENT RNA28S 35 218750 RNA28S||ENST2 GENE2 diff --git a/tests/test_cwl_yaml_adapter.py b/tests/test_cwl_yaml_adapter.py new file mode 100644 index 0000000..6d861a4 --- /dev/null +++ b/tests/test_cwl_yaml_adapter.py @@ -0,0 +1,35 @@ +from __future__ import annotations + +from pathlib import Path + +from workflow.config_adapter import map_cwl_job_to_config, merge_cwl_job_yaml_into_config + + +def test_map_cwl_job_to_config_normalizes_path_objects() -> None: + base = Path('tests/fixtures/mini') + mapped = map_cwl_job_to_config( + { + 'barcode1r1FastqGz': {'class': 'File', 'path': 'source/ip.preRmDup.sam.mini.gz'}, + 'bowtie2_db': {'class': 'Directory', 'path': 'refs/bowtie2_index'}, + 'se_or_pe': 'SE', + 'unrelated': 123, + }, + base_dir=base, + ) + + assert mapped['barcode1r1FastqGz'].endswith('tests/fixtures/mini/source/ip.preRmDup.sam.mini.gz') + assert mapped['bowtie2_db'].endswith('tests/fixtures/mini/refs/bowtie2_index') + assert mapped['se_or_pe'] == 'SE' + assert 'unrelated' not in mapped + + +def test_merge_cwl_job_yaml_into_config_uses_fixture() -> None: + merged = merge_cwl_job_yaml_into_config( + {'dataset': 'base_dataset', 'mini_validation': {'source_sam_gz': 'x'}}, + Path('tests/fixtures/mini/cwl_job_example.yaml'), + ) + + assert merged['dataset'] == 'fixture_dataset' + assert merged['barcode1r1FastqGz'].endswith('tests/fixtures/mini/source/ip.preRmDup.sam.mini.gz') + assert merged['barcode1rmRepBam'] == '/abs/path/ip.bam' + assert merged['mini_validation']['source_sam_gz'] == 'x' diff --git a/tests/test_mini_fixtures_manifest.py b/tests/test_mini_fixtures_manifest.py new file mode 100644 index 0000000..78c77cd --- /dev/null +++ b/tests/test_mini_fixtures_manifest.py @@ -0,0 +1,33 @@ +from __future__ import annotations + +import csv +import hashlib +from pathlib import Path + + +def sha256_file(path: Path) -> str: + h = hashlib.sha256() + with path.open('rb') as fh: + while True: + chunk = fh.read(1024 * 1024) + if not chunk: + break + h.update(chunk) + return h.hexdigest() + + +def test_mini_fixture_manifest_matches_files() -> None: + manifest = Path('tests/fixtures/mini/manifest.tsv') + assert manifest.exists(), 'manifest missing; run scripts/make_mini_fixtures.py' + + with manifest.open('r', encoding='utf-8') as fh: + reader = csv.DictReader(fh, delimiter='\t') + rows = list(reader) + + assert rows, 'manifest contains no fixture rows' + + for row in rows: + path = Path(row['relative_path']) + assert path.exists(), f'missing fixture file: {path}' + assert str(path.stat().st_size) == row['size_bytes'], f'size mismatch: {path}' + assert sha256_file(path) == row['sha256'], f'sha256 mismatch: {path}' diff --git a/tests/test_parse_se_python.py b/tests/test_parse_se_python.py new file mode 100644 index 0000000..dcb60fb --- /dev/null +++ b/tests/test_parse_se_python.py @@ -0,0 +1,49 @@ +from __future__ import annotations + +import subprocess +from pathlib import Path + + +REPO = Path(__file__).resolve().parents[1] +PARSER = REPO / 'bin/python/parse_bowtie2_output_realtime_includemultifamily_SE.py' + + +def test_parse_se_runs_on_sam_input(tmp_path: Path) -> None: + filelist = tmp_path / 'filelist.tsv' + filelist.write_text( + 'ENST0001|ENST0002\tENSGX\tGENEX\tRNA28S\tx\n' + 'ENST0003\tENSGY\tGENEY\tRNA5S\tx\n', + encoding='utf-8', + ) + + sam = tmp_path / 'in.sam' + sam.write_text( + '@HD\tVN:1.0\n' + 'READ1_AA\t0\tENST0001\t1\t255\t10M\t*\t0\t0\tAAAAAAAAAA\tIIIIIIIIII\tAS:i:5\n' + 'READ1_AA\t0\tENST0002\t1\t255\t10M\t*\t0\t0\tAAAAAAAAAA\tIIIIIIIIII\tAS:i:5\n' + 'READ2_TT\t16\tENST0003\t1\t255\t10M\t*\t0\t0\tCCCCCCCCCC\tIIIIIIIIII\tAS:i:4\n', + encoding='utf-8', + ) + + out = tmp_path / 'out.sam' + subprocess.run( + [ + 'python3', + str(PARSER), + 'dummy.fastq.gz', + 'dummy_bowtie_idx', + str(out), + str(filelist), + '--sam-input', + str(sam), + ], + check=True, + ) + + assert out.exists() + assert (tmp_path / 'out.sam.done').exists() + assert (tmp_path / 'out.sam.multimapping_deleted').exists() + + body = [l for l in out.read_text(encoding='utf-8').splitlines() if not l.startswith('@')] + assert body, 'expected emitted mapped lines' + assert all('ZZ:Z:' in l for l in body) diff --git a/tests/test_split_merge_python_expected.py b/tests/test_split_merge_python_expected.py new file mode 100644 index 0000000..455aa2a --- /dev/null +++ b/tests/test_split_merge_python_expected.py @@ -0,0 +1,60 @@ +from __future__ import annotations + +import csv +import gzip +import hashlib +import subprocess +from pathlib import Path + + +REPO = Path(__file__).resolve().parents[1] +PY_SPLIT = REPO / 'bin/python/split_bam_to_subfiles_SEorPE.py' +PY_MERGE = REPO / 'bin/python/merge_multiple_parsed_files.simplified_20191022.py' + + +def sha256_file(path: Path) -> str: + h = hashlib.sha256() + with path.open('rb') as fh: + while True: + chunk = fh.read(1024 * 1024) + if not chunk: + break + h.update(chunk) + return h.hexdigest() + + +def unzip_to(src_gz: Path, dst: Path) -> None: + with gzip.open(src_gz, 'rt', encoding='utf-8', errors='replace') as src, dst.open('w', encoding='utf-8') as out: + out.write(src.read()) + + +def test_split_matches_expected_manifest(tmp_path: Path) -> None: + sam_src = REPO / 'tests/fixtures/mini/source/ip.preRmDup.sam.mini.gz' + sam = tmp_path / 'ip.preRmDup.sam.mini.sam' + unzip_to(sam_src, sam) + + out_dir = tmp_path / 'split' + out_dir.mkdir() + subprocess.run(['python3', str(PY_SPLIT), str(sam), 'SE'], cwd=out_dir, check=True) + + manifest = REPO / 'tests/fixtures/mini/expected/split.expected.manifest.tsv' + with manifest.open('r', encoding='utf-8') as fh: + rows = list(csv.DictReader(fh, delimiter='\t')) + + for row in rows: + p = out_dir / row['filename'] + assert p.exists(), f'missing file: {p.name}' + assert str(p.stat().st_size) == row['size_bytes'] + assert sha256_file(p) == row['sha256'] + + +def test_merge_matches_expected_output(tmp_path: Path) -> None: + in1 = REPO / 'tests/fixtures/mini/source/merge_input_1.parsed_v2.txt' + in2 = REPO / 'tests/fixtures/mini/source/merge_input_2.parsed_v2.txt' + expected = REPO / 'tests/fixtures/mini/expected/merged.expected.parsed' + produced = tmp_path / 'merged.parsed' + + subprocess.run(['python3', str(PY_MERGE), str(produced), str(in1), str(in2)], check=True) + + assert produced.stat().st_size == expected.stat().st_size + assert sha256_file(produced) == sha256_file(expected) diff --git a/wf/eCLIP_repelement_PE b/wf/eCLIP_repelement_PE deleted file mode 100755 index eff6064..0000000 --- a/wf/eCLIP_repelement_PE +++ /dev/null @@ -1,275 +0,0 @@ -#!/usr/bin/env bash - - -#set -o xtrace - - -# # Debugging settings -# #################### -# Print commands and their arguments as they are executed. -#set -o xtrace -#Exit immediately if a command exits with a non-zero status. -set -o errexit -#the return value of a pipeline is the status of -#the last command to exit with a non-zero status, -#or zero if no command exited with a non-zero status -set -o pipefail - - - -CWLRUNNER=cwltool # either: cwltool cwltoil or cwltorq -WORKFLOW=wf_ecliprepmap_pe # the prefix of the workflow cwl (without .cwl extension) - -PIPELINE=REPELEMENTMAPPING # workflow prefix -PIPELINEVERSION=1.0.0 # workflow version - -FULLPATHCOMMAND=$0 -JOB_FILE=$1 - -# get the ecliprepmap yaml file name without extension -JOB_BASENAME=${JOB_FILE##*/} -JOB_NAME="${JOB_BASENAME%.*}" -# TODO the same could be achieved with: JOB_NAME=$(basename "${JOB_FILE}" ".rnaseq") - -# if no parameter, initialize current directory with template files and show help -################# -# FIXME following this usage for * instead of *.cwl does not work -if [[ -z ${JOB_FILE} ]]; then - echo "welcome to rnaseq (torque version)" - # - wf_rnaseqse.usage - exit 0 -fi - - -# from now on, parameter was given (jobfile): -############################################# - -# get the cwlrunner and workflow type -##################################### - -# TODO: harcoded in CWLRUNNER and WORKFLOW vars above - -# remove the path at the beginning of full path command -# CWLRUNNER_WORKFLOW=${FULLPATHCOMMAND##*/} -# echo "FULL PATH COMMAND" -# echo $CWLRUNNER_WORKFLOW; - -# if command was 'rnaseqse' then will use the default cwlrunner and the default workflow -# if [ $CWLRUNNER_WORKFLOW = 'rnaseqse' ] -# then -# #echo default cwlrunner and workflow -# CWLRUNNER=cwltool -# WORKFLOW=wf_rnaseqcore_se -# else -# # remove cwlto+6 chars from beginning -# # WORKFLOW=${CWLRUNNER_WORKFLOW#cwlto??} # ?? matches either OL or IL or RQ # TODO caution: no .cwl extension ! -# WORKFLOW=wf_rnaseqcore_se_cwltoil -# # get those cwl+6 chars that were removed above -# # CWLRUNNER=${CWLRUNNER_WORKFLOW%$WORKFLOW} # FIXME: don't know what this does -# CWLRUNNER=cwltorq # FIXME: hardcoded for now... -# if [ -z $WORKFLOW ] -# then -# WORKFLOW=wf_rnaseqcore_se_cwltorq -# fi -# fi - - - -# initialize directories -######################## - -# get this script full directory name no matter where it is called from -WORKFLOW_HOME="$( cd "$( dirname "${BASH_SOURCE[0]}" )/../" && pwd )" # TODO move this one level up (include final /../) - - -#export ECLIPREPMAP_DOCUMENTATION=${WORKFLOW_HOME}/../documentation -#export WORKFLOW_HOME=${WORKFLOW_HOME}/../init -# get the eclipjob.yaml full directory name -JOB_HOME="$( cd "$( dirname "${BASH_SOURCE[1]}" )" && pwd )" - -# create temporary intermediates directory # TODO only needed for cwltool -INTERM=${JOB_HOME}/${JOB_NAME}/.tmp/cwltool_interm -mkdir -p ${INTERM} - -# create temp work directory # TODO should be absolute path for toil, is it ? -WORKDIR=${JOB_HOME}/${JOB_NAME}/.tmp/workdir -mkdir -p ${WORKDIR} -echo "WORKDIR: ${WORKDIR}" - -# define jobstore directory, but do not create it -JOBSTORE=${JOB_HOME}/${JOB_NAME}/.tmp/cwltoil_jobstore -#mkdir -p ${JOBSTORE}; -JOBSTOREEXITSLETSRESTART= -# if jobstore exists set arguments to restart toil -if [ -d $JOBSTORE ] -then echo jobstore exists; JOBSTOREEXITSLETSRESTART="--restart" -fi - -# create outdir directory -OUTDIR=${JOB_HOME}/${JOB_NAME}/.tmp/outdir -mkdir -p ${OUTDIR} - -# create toillogs directory -TOILLOGS=${JOB_HOME}/${JOB_NAME}/.tmp/toillogs -mkdir -p ${TOILLOGS} - -# create results directory -RESULTSDIR=${JOB_HOME}/${JOB_NAME}/results -mkdir -p ${RESULTSDIR} - - -# copy job file and log ecliprepmap version in output directory -######################################################### -# TODO softcode version nb -JOBCOPYFILEPATH="${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_INPUT.yaml" -cp ${JOB_FILE} ${JOBCOPYFILEPATH} -chmod -x ${JOBCOPYFILEPATH} -touch ${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_VERSION-${PIPELINEVERSION} -# TODO: does not do anything -# touch ${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_PIPELINE-rnaseqse -touch ${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_WORKFLOW-${WORKFLOW} - -echo -echo ======================================================================== -echo PATHS -echo ======================================================================== -echo -echo "PATH:" $PATH -echo "which conda:" `which conda` -# echo "which rnaseqse:" `which rnaseqse` -echo "JOB_HOME:" $JOB_HOME -echo "OUTDIR:" $OUTDIR -echo "RESULTSDIR:" $RESULTSDIR - - -# execute -########## -echo -echo ======================================================================== -echo EXECUTING CWL JOB -echo ======================================================================== -echo - -if [ $CWLRUNNER == cwltool ] -then - cwltool --debug \ - --outdir ${OUTDIR} \ - --cachedir ${INTERM} \ - ${WORKFLOW_HOME}/cwl/${WORKFLOW}.cwl \ - ${JOB_FILE} \ - 2>${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_LOG.txt | tee ${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_OUTPUT.json - - #2>&1 | tee ${JOB_HOME}/${JOB_NAME}/ECLIPREPMAP-JOB-LOG.txt - -elif [ $CWLRUNNER == cwltoil -o $CWLRUNNER == cwltorq ] -then - if [ $CWLRUNNER == cwltorq ]; then - BATCHSYSTEM="torque --disableCaching" - export TOIL_TORQUE_ARGS="-q home-yeo" # home-yeo - export TOIL_TORQUE_REQS="walltime=12:00:00" # walltime=8:00:00 - else - BATCHSYSTEM="single_machine" - fi - echo - echo CWLRUNNER : $CWLRUNNER - echo BATCHSYSTEM: $BATCHSYSTEM - echo - # For more robust logging options (default is INFO), - # use --logDebug - # or more generally, use --logLevel=, which may be set to either - # OFF (or CRITICAL), ERROR, WARN (or WARNING), INFO or DEBUG. - # Logs can be directed to a file with --logFile=. - # --logLevel DEBUG \ - # --logLevel=critical \ - # --logLevel=error \ - # --logLevel=warning \ - # --logLevel=info \ - # --logLevel=debug \ - # - # The defaults for cores (1), disk (2G), and memory (2G), can all be changed - # using --defaultCores, --defaultDisk, and --defaultMemory. - # Standard suffixes like K, Ki, M, Mi, G or Gi are supported. - - -# ${WORKFLOW_HOME}/repo/cwl/${WORKFLOW}.cwl \ - - cwltoil ${JOBSTOREEXITSLETSRESTART} \ - --batchSystem ${BATCHSYSTEM} \ - --stats \ - --writeLogs ${TOILLOGS} \ - --jobStore file:${JOBSTORE} \ - --workDir ${WORKDIR} \ - --outdir ${OUTDIR} \ - --clean never \ - --cleanWorkDir onSuccess \ - --logInfo \ - --clusterStats \ - --realTimeLogging \ - --stats \ - --defaultMemory 64.0G \ - --no-container \ - ${WORKFLOW_HOME}/cwl/${WORKFLOW}.cwl \ - ${JOB_FILE} \ - 2>${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_LOG.txt | tee ${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_OUTPUT.json -fi - -#--defaultCores 1.0 \ -#--defaultMemory 8.0G \ -#--clusterStats \ -#--logDebug \ -#--realTimeLogging \ -#--stats \ -#--maxCores= --maxMemory=8.0Ei --maxDisk 8.0Ei \ -#--maxLogFileSize -1 \ - -# error when using option --maxLogFileSize -1 -# h2b = lambda x : human2bytes(str(x)) -#File "/projects/ps-yeolab/software/eclipconda/envs/ecliplesspeakstoil/lib/python2.7/site-packages/bd2k/util/humanize.py", line 111, in human2bytes -# num = float(num) - -############################################ -#rmdir ${JOB_HOME}/${JOB_NAME}/intermediates?????? - - -# Stop printing commands and their arguments as they are executed. -set +x - - - -echo -echo ======================================================================== -echo UPDATING RESULTS -echo ======================================================================== -echo -echo tmp output dir: ${OUTDIR} -echo results dir: ${RESULTSDIR} - -mv ${OUTDIR}/* ${RESULTSDIR}/ - -echo "Moving tmp to outdir done." - -cd ${JOB_HOME}/${JOB_NAME} - -if grep -q "permanentFail" ${PIPELINE}_${JOB_NAME}_LOG.txt -then - echo "Not good (failures found, check log). Keeping intermediates for now." # code if found -else - echo "All good (no failures)! Removing intermediate tmp folder." # code if not found - rm -rf ./.tmp/; - rm -rf ${RESULTSDIR}/out_tmp* ${RESULTSDIR}/tmp* -fi - -# TODO this now commented out, so to finish inside job folder -#cd - - -# Stop printing commands and their arguments as they are executed. - -echo -echo "*********************************************************************" -echo "* SUCCESS: ${JOB_NAME} PROCESSED WITH ${PIPELINEVERSION} *" -echo "*********************************************************************" -echo -echo Bye! -echo - diff --git a/wf/eCLIP_repelement_PE_singleNode b/wf/eCLIP_repelement_PE_singleNode deleted file mode 100755 index 0d6de73..0000000 --- a/wf/eCLIP_repelement_PE_singleNode +++ /dev/null @@ -1,275 +0,0 @@ -#!/usr/bin/env bash - - -#set -o xtrace - - -# # Debugging settings -# #################### -# Print commands and their arguments as they are executed. -#set -o xtrace -#Exit immediately if a command exits with a non-zero status. -set -o errexit -#the return value of a pipeline is the status of -#the last command to exit with a non-zero status, -#or zero if no command exited with a non-zero status -set -o pipefail - - - -CWLRUNNER=cwltool # either: cwltool cwltoil or cwltorq -WORKFLOW=wf_ecliprepmap_pe # the prefix of the workflow cwl (without .cwl extension) - -PIPELINE=REPELEMENTMAPPING # workflow prefix -PIPELINEVERSION=0.1.0 # workflow version - -FULLPATHCOMMAND=$0 -JOB_FILE=$1 - -# get the ecliprepmap yaml file name without extension -JOB_BASENAME=${JOB_FILE##*/} -JOB_NAME="${JOB_BASENAME%.*}" -# TODO the same could be achieved with: JOB_NAME=$(basename "${JOB_FILE}" ".rnaseq") - -# if no parameter, initialize current directory with template files and show help -################# -# FIXME following this usage for * instead of *.cwl does not work -if [[ -z ${JOB_FILE} ]]; then - echo "welcome to rnaseq (torque version)" - # - wf_rnaseqse.usage - exit 0 -fi - - -# from now on, parameter was given (jobfile): -############################################# - -# get the cwlrunner and workflow type -##################################### - -# TODO: harcoded in CWLRUNNER and WORKFLOW vars above - -# remove the path at the beginning of full path command -# CWLRUNNER_WORKFLOW=${FULLPATHCOMMAND##*/} -# echo "FULL PATH COMMAND" -# echo $CWLRUNNER_WORKFLOW; - -# if command was 'rnaseqse' then will use the default cwlrunner and the default workflow -# if [ $CWLRUNNER_WORKFLOW = 'rnaseqse' ] -# then -# #echo default cwlrunner and workflow -# CWLRUNNER=cwltool -# WORKFLOW=wf_rnaseqcore_se -# else -# # remove cwlto+6 chars from beginning -# # WORKFLOW=${CWLRUNNER_WORKFLOW#cwlto??} # ?? matches either OL or IL or RQ # TODO caution: no .cwl extension ! -# WORKFLOW=wf_rnaseqcore_se_cwltoil -# # get those cwl+6 chars that were removed above -# # CWLRUNNER=${CWLRUNNER_WORKFLOW%$WORKFLOW} # FIXME: don't know what this does -# CWLRUNNER=cwltorq # FIXME: hardcoded for now... -# if [ -z $WORKFLOW ] -# then -# WORKFLOW=wf_rnaseqcore_se_cwltorq -# fi -# fi - - - -# initialize directories -######################## - -# get this script full directory name no matter where it is called from -WORKFLOW_HOME="$( cd "$( dirname "${BASH_SOURCE[0]}" )/../" && pwd )" # TODO move this one level up (include final /../) - - -#export ECLIPREPMAP_DOCUMENTATION=${WORKFLOW_HOME}/../documentation -#export WORKFLOW_HOME=${WORKFLOW_HOME}/../init -# get the eclipjob.yaml full directory name -JOB_HOME="$( cd "$( dirname "${BASH_SOURCE[1]}" )" && pwd )" - -# create temporary intermediates directory # TODO only needed for cwltool -INTERM=${JOB_HOME}/${JOB_NAME}/.tmp/cwltool_interm -mkdir -p ${INTERM} - -# create temp work directory # TODO should be absolute path for toil, is it ? -WORKDIR=${JOB_HOME}/${JOB_NAME}/.tmp/workdir -mkdir -p ${WORKDIR} -echo "WORKDIR: ${WORKDIR}" - -# define jobstore directory, but do not create it -JOBSTORE=${JOB_HOME}/${JOB_NAME}/.tmp/cwltoil_jobstore -#mkdir -p ${JOBSTORE}; -JOBSTOREEXITSLETSRESTART= -# if jobstore exists set arguments to restart toil -if [ -d $JOBSTORE ] -then echo jobstore exists; JOBSTOREEXITSLETSRESTART="--restart" -fi - -# create outdir directory -OUTDIR=${JOB_HOME}/${JOB_NAME}/.tmp/outdir -mkdir -p ${OUTDIR} - -# create toillogs directory -TOILLOGS=${JOB_HOME}/${JOB_NAME}/.tmp/toillogs -mkdir -p ${TOILLOGS} - -# create results directory -RESULTSDIR=${JOB_HOME}/${JOB_NAME}/results -mkdir -p ${RESULTSDIR} - - -# copy job file and log ecliprepmap version in output directory -######################################################### -# TODO softcode version nb -JOBCOPYFILEPATH="${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_INPUT.yaml" -cp ${JOB_FILE} ${JOBCOPYFILEPATH} -chmod -x ${JOBCOPYFILEPATH} -touch ${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_VERSION-${PIPELINEVERSION} -# TODO: does not do anything -# touch ${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_PIPELINE-rnaseqse -touch ${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_WORKFLOW-${WORKFLOW} - -echo -echo ======================================================================== -echo PATHS -echo ======================================================================== -echo -echo "PATH:" $PATH -echo "which conda:" `which conda` -# echo "which rnaseqse:" `which rnaseqse` -echo "JOB_HOME:" $JOB_HOME -echo "OUTDIR:" $OUTDIR -echo "RESULTSDIR:" $RESULTSDIR - - -# execute -########## -echo -echo ======================================================================== -echo EXECUTING CWL JOB -echo ======================================================================== -echo - -if [ $CWLRUNNER == cwltool ] -then - cwltool --debug \ - --outdir ${OUTDIR} \ - --cachedir ${INTERM} \ - ${WORKFLOW_HOME}/cwl/${WORKFLOW}.cwl \ - ${JOB_FILE} \ - 2>${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_LOG.txt | tee ${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_OUTPUT.json - - #2>&1 | tee ${JOB_HOME}/${JOB_NAME}/ECLIPREPMAP-JOB-LOG.txt - -elif [ $CWLRUNNER == cwltoil -o $CWLRUNNER == cwltorq ] -then - if [ $CWLRUNNER == cwltorq ]; then - BATCHSYSTEM="torque --disableCaching" - export TOIL_TORQUE_ARGS="-q home-yeo" # home-yeo - export TOIL_TORQUE_REQS="walltime=12:00:00" # walltime=8:00:00 - else - BATCHSYSTEM="single_machine" - fi - echo - echo CWLRUNNER : $CWLRUNNER - echo BATCHSYSTEM: $BATCHSYSTEM - echo - # For more robust logging options (default is INFO), - # use --logDebug - # or more generally, use --logLevel=, which may be set to either - # OFF (or CRITICAL), ERROR, WARN (or WARNING), INFO or DEBUG. - # Logs can be directed to a file with --logFile=. - # --logLevel DEBUG \ - # --logLevel=critical \ - # --logLevel=error \ - # --logLevel=warning \ - # --logLevel=info \ - # --logLevel=debug \ - # - # The defaults for cores (1), disk (2G), and memory (2G), can all be changed - # using --defaultCores, --defaultDisk, and --defaultMemory. - # Standard suffixes like K, Ki, M, Mi, G or Gi are supported. - - -# ${WORKFLOW_HOME}/repo/cwl/${WORKFLOW}.cwl \ - - cwltoil ${JOBSTOREEXITSLETSRESTART} \ - --batchSystem ${BATCHSYSTEM} \ - --stats \ - --writeLogs ${TOILLOGS} \ - --jobStore file:${JOBSTORE} \ - --workDir ${WORKDIR} \ - --outdir ${OUTDIR} \ - --clean never \ - --cleanWorkDir onSuccess \ - --logInfo \ - --clusterStats \ - --realTimeLogging \ - --stats \ - --defaultMemory 64.0G \ - --no-container \ - ${WORKFLOW_HOME}/cwl/${WORKFLOW}.cwl \ - ${JOB_FILE} \ - 2>${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_LOG.txt | tee ${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_OUTPUT.json -fi - -#--defaultCores 1.0 \ -#--defaultMemory 8.0G \ -#--clusterStats \ -#--logDebug \ -#--realTimeLogging \ -#--stats \ -#--maxCores= --maxMemory=8.0Ei --maxDisk 8.0Ei \ -#--maxLogFileSize -1 \ - -# error when using option --maxLogFileSize -1 -# h2b = lambda x : human2bytes(str(x)) -#File "/projects/ps-yeolab/software/eclipconda/envs/ecliplesspeakstoil/lib/python2.7/site-packages/bd2k/util/humanize.py", line 111, in human2bytes -# num = float(num) - -############################################ -#rmdir ${JOB_HOME}/${JOB_NAME}/intermediates?????? - - -# Stop printing commands and their arguments as they are executed. -set +x - - - -echo -echo ======================================================================== -echo UPDATING RESULTS -echo ======================================================================== -echo -echo tmp output dir: ${OUTDIR} -echo results dir: ${RESULTSDIR} - -mv ${OUTDIR}/* ${RESULTSDIR}/ - -echo "Moving tmp to outdir done." - -cd ${JOB_HOME}/${JOB_NAME} - -if grep -q "permanentFail" ${PIPELINE}_${JOB_NAME}_LOG.txt -then - echo "Not good (failures found, check log). Keeping intermediates for now." # code if found -else - echo "All good (no failures)! Removing intermediate tmp folder." # code if not found - rm -rf ./.tmp/; - rm -rf ${RESULTSDIR}/out_tmp* ${RESULTSDIR}/tmp* -fi - -# TODO this now commented out, so to finish inside job folder -#cd - - -# Stop printing commands and their arguments as they are executed. - -echo -echo "*********************************************************************" -echo "* SUCCESS: ${JOB_NAME} PROCESSED WITH ${PIPELINEVERSION} *" -echo "*********************************************************************" -echo -echo Bye! -echo - diff --git a/wf/eCLIP_repelement_SE b/wf/eCLIP_repelement_SE deleted file mode 100755 index 9a467ea..0000000 --- a/wf/eCLIP_repelement_SE +++ /dev/null @@ -1,275 +0,0 @@ -#!/usr/bin/env bash - - -#set -o xtrace - - -# # Debugging settings -# #################### -# Print commands and their arguments as they are executed. -#set -o xtrace -#Exit immediately if a command exits with a non-zero status. -set -o errexit -#the return value of a pipeline is the status of -#the last command to exit with a non-zero status, -#or zero if no command exited with a non-zero status -set -o pipefail - - - -CWLRUNNER=cwltool # either: cwltool cwltoil or cwltorq -WORKFLOW=wf_ecliprepmap_se # the prefix of the workflow cwl (without .cwl extension) - -PIPELINE=REPELEMENTMAPPING # workflow prefix -PIPELINEVERSION=1.0.0 # workflow version - -FULLPATHCOMMAND=$0 -JOB_FILE=$1 - -# get the ecliprepmap yaml file name without extension -JOB_BASENAME=${JOB_FILE##*/} -JOB_NAME="${JOB_BASENAME%.*}" -# TODO the same could be achieved with: JOB_NAME=$(basename "${JOB_FILE}" ".rnaseq") - -# if no parameter, initialize current directory with template files and show help -################# -# FIXME following this usage for * instead of *.cwl does not work -if [[ -z ${JOB_FILE} ]]; then - echo "welcome to rnaseq (torque version)" - # - wf_rnaseqse.usage - exit 0 -fi - - -# from now on, parameter was given (jobfile): -############################################# - -# get the cwlrunner and workflow type -##################################### - -# TODO: harcoded in CWLRUNNER and WORKFLOW vars above - -# remove the path at the beginning of full path command -# CWLRUNNER_WORKFLOW=${FULLPATHCOMMAND##*/} -# echo "FULL PATH COMMAND" -# echo $CWLRUNNER_WORKFLOW; - -# if command was 'rnaseqse' then will use the default cwlrunner and the default workflow -# if [ $CWLRUNNER_WORKFLOW = 'rnaseqse' ] -# then -# #echo default cwlrunner and workflow -# CWLRUNNER=cwltool -# WORKFLOW=wf_rnaseqcore_se -# else -# # remove cwlto+6 chars from beginning -# # WORKFLOW=${CWLRUNNER_WORKFLOW#cwlto??} # ?? matches either OL or IL or RQ # TODO caution: no .cwl extension ! -# WORKFLOW=wf_rnaseqcore_se_cwltoil -# # get those cwl+6 chars that were removed above -# # CWLRUNNER=${CWLRUNNER_WORKFLOW%$WORKFLOW} # FIXME: don't know what this does -# CWLRUNNER=cwltorq # FIXME: hardcoded for now... -# if [ -z $WORKFLOW ] -# then -# WORKFLOW=wf_rnaseqcore_se_cwltorq -# fi -# fi - - - -# initialize directories -######################## - -# get this script full directory name no matter where it is called from -WORKFLOW_HOME="$( cd "$( dirname "${BASH_SOURCE[0]}" )/../" && pwd )" # TODO move this one level up (include final /../) - - -#export ECLIPREPMAP_DOCUMENTATION=${WORKFLOW_HOME}/../documentation -#export WORKFLOW_HOME=${WORKFLOW_HOME}/../init -# get the eclipjob.yaml full directory name -JOB_HOME="$( cd "$( dirname "${BASH_SOURCE[1]}" )" && pwd )" - -# create temporary intermediates directory # TODO only needed for cwltool -INTERM=${JOB_HOME}/${JOB_NAME}/.tmp/cwltool_interm -mkdir -p ${INTERM} - -# create temp work directory # TODO should be absolute path for toil, is it ? -WORKDIR=${JOB_HOME}/${JOB_NAME}/.tmp/workdir -mkdir -p ${WORKDIR} -echo "WORKDIR: ${WORKDIR}" - -# define jobstore directory, but do not create it -JOBSTORE=${JOB_HOME}/${JOB_NAME}/.tmp/cwltoil_jobstore -#mkdir -p ${JOBSTORE}; -JOBSTOREEXITSLETSRESTART= -# if jobstore exists set arguments to restart toil -if [ -d $JOBSTORE ] -then echo jobstore exists; JOBSTOREEXITSLETSRESTART="--restart" -fi - -# create outdir directory -OUTDIR=${JOB_HOME}/${JOB_NAME}/.tmp/outdir -mkdir -p ${OUTDIR} - -# create toillogs directory -TOILLOGS=${JOB_HOME}/${JOB_NAME}/.tmp/toillogs -mkdir -p ${TOILLOGS} - -# create results directory -RESULTSDIR=${JOB_HOME}/${JOB_NAME}/results -mkdir -p ${RESULTSDIR} - - -# copy job file and log ecliprepmap version in output directory -######################################################### -# TODO softcode version nb -JOBCOPYFILEPATH="${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_INPUT.yaml" -cp ${JOB_FILE} ${JOBCOPYFILEPATH} -chmod -x ${JOBCOPYFILEPATH} -touch ${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_VERSION-${PIPELINEVERSION} -# TODO: does not do anything -# touch ${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_PIPELINE-rnaseqse -touch ${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_WORKFLOW-${WORKFLOW} - -echo -echo ======================================================================== -echo PATHS -echo ======================================================================== -echo -echo "PATH:" $PATH -echo "which conda:" `which conda` -# echo "which rnaseqse:" `which rnaseqse` -echo "JOB_HOME:" $JOB_HOME -echo "OUTDIR:" $OUTDIR -echo "RESULTSDIR:" $RESULTSDIR - - -# execute -########## -echo -echo ======================================================================== -echo EXECUTING CWL JOB -echo ======================================================================== -echo - -if [ $CWLRUNNER == cwltool ] -then - cwltool --debug \ - --outdir ${OUTDIR} \ - --cachedir ${INTERM} \ - ${WORKFLOW_HOME}/cwl/${WORKFLOW}.cwl \ - ${JOB_FILE} \ - 2>${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_LOG.txt | tee ${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_OUTPUT.json - - #2>&1 | tee ${JOB_HOME}/${JOB_NAME}/ECLIPREPMAP-JOB-LOG.txt - -elif [ $CWLRUNNER == cwltoil -o $CWLRUNNER == cwltorq ] -then - if [ $CWLRUNNER == cwltorq ]; then - BATCHSYSTEM="torque --disableCaching" - export TOIL_TORQUE_ARGS="-q home-yeo" # home-yeo - export TOIL_TORQUE_REQS="walltime=12:00:00" # walltime=8:00:00 - else - BATCHSYSTEM="single_machine" - fi - echo - echo CWLRUNNER : $CWLRUNNER - echo BATCHSYSTEM: $BATCHSYSTEM - echo - # For more robust logging options (default is INFO), - # use --logDebug - # or more generally, use --logLevel=, which may be set to either - # OFF (or CRITICAL), ERROR, WARN (or WARNING), INFO or DEBUG. - # Logs can be directed to a file with --logFile=. - # --logLevel DEBUG \ - # --logLevel=critical \ - # --logLevel=error \ - # --logLevel=warning \ - # --logLevel=info \ - # --logLevel=debug \ - # - # The defaults for cores (1), disk (2G), and memory (2G), can all be changed - # using --defaultCores, --defaultDisk, and --defaultMemory. - # Standard suffixes like K, Ki, M, Mi, G or Gi are supported. - - -# ${WORKFLOW_HOME}/repo/cwl/${WORKFLOW}.cwl \ - - cwltoil ${JOBSTOREEXITSLETSRESTART} \ - --batchSystem ${BATCHSYSTEM} \ - --stats \ - --writeLogs ${TOILLOGS} \ - --jobStore file:${JOBSTORE} \ - --workDir ${WORKDIR} \ - --outdir ${OUTDIR} \ - --clean never \ - --cleanWorkDir onSuccess \ - --logInfo \ - --clusterStats \ - --realTimeLogging \ - --stats \ - --defaultMemory 64.0G \ - --no-container \ - ${WORKFLOW_HOME}/cwl/${WORKFLOW}.cwl \ - ${JOB_FILE} \ - 2>${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_LOG.txt | tee ${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_OUTPUT.json -fi -# --defaultMemory 64.0G \ -#--defaultCores 1.0 \ -#--defaultMemory 8.0G \ -#--clusterStats \ -#--logDebug \ -#--realTimeLogging \ -#--stats \ -#--maxCores= --maxMemory=8.0Ei --maxDisk 8.0Ei \ -#--maxLogFileSize -1 \ - -# error when using option --maxLogFileSize -1 -# h2b = lambda x : human2bytes(str(x)) -#File "/projects/ps-yeolab/software/eclipconda/envs/ecliplesspeakstoil/lib/python2.7/site-packages/bd2k/util/humanize.py", line 111, in human2bytes -# num = float(num) - -############################################ -#rmdir ${JOB_HOME}/${JOB_NAME}/intermediates?????? - - -# Stop printing commands and their arguments as they are executed. -set +x - - - -echo -echo ======================================================================== -echo UPDATING RESULTS -echo ======================================================================== -echo -echo tmp output dir: ${OUTDIR} -echo results dir: ${RESULTSDIR} - -mv ${OUTDIR}/* ${RESULTSDIR}/ - -echo "Moving tmp to outdir done." - -cd ${JOB_HOME}/${JOB_NAME} - -if grep -q "permanentFail" ${PIPELINE}_${JOB_NAME}_LOG.txt -then - echo "Not good (failures found, check log). Keeping intermediates for now." # code if found -else - echo "All good (no failures)! Removing intermediate tmp folder." # code if not found - rm -rf ./.tmp/; - rm -rf ${RESULTSDIR}/out_tmp* ${RESULTSDIR}/tmp* -fi - -# TODO this now commented out, so to finish inside job folder -#cd - - -# Stop printing commands and their arguments as they are executed. - -echo -echo "*********************************************************************" -echo "* SUCCESS: ${JOB_NAME} PROCESSED WITH ${PIPELINEVERSION} *" -echo "*********************************************************************" -echo -echo Bye! -echo - diff --git a/wf/eCLIP_repelement_SE_singleNode b/wf/eCLIP_repelement_SE_singleNode deleted file mode 100755 index 9a467ea..0000000 --- a/wf/eCLIP_repelement_SE_singleNode +++ /dev/null @@ -1,275 +0,0 @@ -#!/usr/bin/env bash - - -#set -o xtrace - - -# # Debugging settings -# #################### -# Print commands and their arguments as they are executed. -#set -o xtrace -#Exit immediately if a command exits with a non-zero status. -set -o errexit -#the return value of a pipeline is the status of -#the last command to exit with a non-zero status, -#or zero if no command exited with a non-zero status -set -o pipefail - - - -CWLRUNNER=cwltool # either: cwltool cwltoil or cwltorq -WORKFLOW=wf_ecliprepmap_se # the prefix of the workflow cwl (without .cwl extension) - -PIPELINE=REPELEMENTMAPPING # workflow prefix -PIPELINEVERSION=1.0.0 # workflow version - -FULLPATHCOMMAND=$0 -JOB_FILE=$1 - -# get the ecliprepmap yaml file name without extension -JOB_BASENAME=${JOB_FILE##*/} -JOB_NAME="${JOB_BASENAME%.*}" -# TODO the same could be achieved with: JOB_NAME=$(basename "${JOB_FILE}" ".rnaseq") - -# if no parameter, initialize current directory with template files and show help -################# -# FIXME following this usage for * instead of *.cwl does not work -if [[ -z ${JOB_FILE} ]]; then - echo "welcome to rnaseq (torque version)" - # - wf_rnaseqse.usage - exit 0 -fi - - -# from now on, parameter was given (jobfile): -############################################# - -# get the cwlrunner and workflow type -##################################### - -# TODO: harcoded in CWLRUNNER and WORKFLOW vars above - -# remove the path at the beginning of full path command -# CWLRUNNER_WORKFLOW=${FULLPATHCOMMAND##*/} -# echo "FULL PATH COMMAND" -# echo $CWLRUNNER_WORKFLOW; - -# if command was 'rnaseqse' then will use the default cwlrunner and the default workflow -# if [ $CWLRUNNER_WORKFLOW = 'rnaseqse' ] -# then -# #echo default cwlrunner and workflow -# CWLRUNNER=cwltool -# WORKFLOW=wf_rnaseqcore_se -# else -# # remove cwlto+6 chars from beginning -# # WORKFLOW=${CWLRUNNER_WORKFLOW#cwlto??} # ?? matches either OL or IL or RQ # TODO caution: no .cwl extension ! -# WORKFLOW=wf_rnaseqcore_se_cwltoil -# # get those cwl+6 chars that were removed above -# # CWLRUNNER=${CWLRUNNER_WORKFLOW%$WORKFLOW} # FIXME: don't know what this does -# CWLRUNNER=cwltorq # FIXME: hardcoded for now... -# if [ -z $WORKFLOW ] -# then -# WORKFLOW=wf_rnaseqcore_se_cwltorq -# fi -# fi - - - -# initialize directories -######################## - -# get this script full directory name no matter where it is called from -WORKFLOW_HOME="$( cd "$( dirname "${BASH_SOURCE[0]}" )/../" && pwd )" # TODO move this one level up (include final /../) - - -#export ECLIPREPMAP_DOCUMENTATION=${WORKFLOW_HOME}/../documentation -#export WORKFLOW_HOME=${WORKFLOW_HOME}/../init -# get the eclipjob.yaml full directory name -JOB_HOME="$( cd "$( dirname "${BASH_SOURCE[1]}" )" && pwd )" - -# create temporary intermediates directory # TODO only needed for cwltool -INTERM=${JOB_HOME}/${JOB_NAME}/.tmp/cwltool_interm -mkdir -p ${INTERM} - -# create temp work directory # TODO should be absolute path for toil, is it ? -WORKDIR=${JOB_HOME}/${JOB_NAME}/.tmp/workdir -mkdir -p ${WORKDIR} -echo "WORKDIR: ${WORKDIR}" - -# define jobstore directory, but do not create it -JOBSTORE=${JOB_HOME}/${JOB_NAME}/.tmp/cwltoil_jobstore -#mkdir -p ${JOBSTORE}; -JOBSTOREEXITSLETSRESTART= -# if jobstore exists set arguments to restart toil -if [ -d $JOBSTORE ] -then echo jobstore exists; JOBSTOREEXITSLETSRESTART="--restart" -fi - -# create outdir directory -OUTDIR=${JOB_HOME}/${JOB_NAME}/.tmp/outdir -mkdir -p ${OUTDIR} - -# create toillogs directory -TOILLOGS=${JOB_HOME}/${JOB_NAME}/.tmp/toillogs -mkdir -p ${TOILLOGS} - -# create results directory -RESULTSDIR=${JOB_HOME}/${JOB_NAME}/results -mkdir -p ${RESULTSDIR} - - -# copy job file and log ecliprepmap version in output directory -######################################################### -# TODO softcode version nb -JOBCOPYFILEPATH="${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_INPUT.yaml" -cp ${JOB_FILE} ${JOBCOPYFILEPATH} -chmod -x ${JOBCOPYFILEPATH} -touch ${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_VERSION-${PIPELINEVERSION} -# TODO: does not do anything -# touch ${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_PIPELINE-rnaseqse -touch ${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_WORKFLOW-${WORKFLOW} - -echo -echo ======================================================================== -echo PATHS -echo ======================================================================== -echo -echo "PATH:" $PATH -echo "which conda:" `which conda` -# echo "which rnaseqse:" `which rnaseqse` -echo "JOB_HOME:" $JOB_HOME -echo "OUTDIR:" $OUTDIR -echo "RESULTSDIR:" $RESULTSDIR - - -# execute -########## -echo -echo ======================================================================== -echo EXECUTING CWL JOB -echo ======================================================================== -echo - -if [ $CWLRUNNER == cwltool ] -then - cwltool --debug \ - --outdir ${OUTDIR} \ - --cachedir ${INTERM} \ - ${WORKFLOW_HOME}/cwl/${WORKFLOW}.cwl \ - ${JOB_FILE} \ - 2>${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_LOG.txt | tee ${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_OUTPUT.json - - #2>&1 | tee ${JOB_HOME}/${JOB_NAME}/ECLIPREPMAP-JOB-LOG.txt - -elif [ $CWLRUNNER == cwltoil -o $CWLRUNNER == cwltorq ] -then - if [ $CWLRUNNER == cwltorq ]; then - BATCHSYSTEM="torque --disableCaching" - export TOIL_TORQUE_ARGS="-q home-yeo" # home-yeo - export TOIL_TORQUE_REQS="walltime=12:00:00" # walltime=8:00:00 - else - BATCHSYSTEM="single_machine" - fi - echo - echo CWLRUNNER : $CWLRUNNER - echo BATCHSYSTEM: $BATCHSYSTEM - echo - # For more robust logging options (default is INFO), - # use --logDebug - # or more generally, use --logLevel=, which may be set to either - # OFF (or CRITICAL), ERROR, WARN (or WARNING), INFO or DEBUG. - # Logs can be directed to a file with --logFile=. - # --logLevel DEBUG \ - # --logLevel=critical \ - # --logLevel=error \ - # --logLevel=warning \ - # --logLevel=info \ - # --logLevel=debug \ - # - # The defaults for cores (1), disk (2G), and memory (2G), can all be changed - # using --defaultCores, --defaultDisk, and --defaultMemory. - # Standard suffixes like K, Ki, M, Mi, G or Gi are supported. - - -# ${WORKFLOW_HOME}/repo/cwl/${WORKFLOW}.cwl \ - - cwltoil ${JOBSTOREEXITSLETSRESTART} \ - --batchSystem ${BATCHSYSTEM} \ - --stats \ - --writeLogs ${TOILLOGS} \ - --jobStore file:${JOBSTORE} \ - --workDir ${WORKDIR} \ - --outdir ${OUTDIR} \ - --clean never \ - --cleanWorkDir onSuccess \ - --logInfo \ - --clusterStats \ - --realTimeLogging \ - --stats \ - --defaultMemory 64.0G \ - --no-container \ - ${WORKFLOW_HOME}/cwl/${WORKFLOW}.cwl \ - ${JOB_FILE} \ - 2>${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_LOG.txt | tee ${JOB_HOME}/${JOB_NAME}/${PIPELINE}_${JOB_NAME}_OUTPUT.json -fi -# --defaultMemory 64.0G \ -#--defaultCores 1.0 \ -#--defaultMemory 8.0G \ -#--clusterStats \ -#--logDebug \ -#--realTimeLogging \ -#--stats \ -#--maxCores= --maxMemory=8.0Ei --maxDisk 8.0Ei \ -#--maxLogFileSize -1 \ - -# error when using option --maxLogFileSize -1 -# h2b = lambda x : human2bytes(str(x)) -#File "/projects/ps-yeolab/software/eclipconda/envs/ecliplesspeakstoil/lib/python2.7/site-packages/bd2k/util/humanize.py", line 111, in human2bytes -# num = float(num) - -############################################ -#rmdir ${JOB_HOME}/${JOB_NAME}/intermediates?????? - - -# Stop printing commands and their arguments as they are executed. -set +x - - - -echo -echo ======================================================================== -echo UPDATING RESULTS -echo ======================================================================== -echo -echo tmp output dir: ${OUTDIR} -echo results dir: ${RESULTSDIR} - -mv ${OUTDIR}/* ${RESULTSDIR}/ - -echo "Moving tmp to outdir done." - -cd ${JOB_HOME}/${JOB_NAME} - -if grep -q "permanentFail" ${PIPELINE}_${JOB_NAME}_LOG.txt -then - echo "Not good (failures found, check log). Keeping intermediates for now." # code if found -else - echo "All good (no failures)! Removing intermediate tmp folder." # code if not found - rm -rf ./.tmp/; - rm -rf ${RESULTSDIR}/out_tmp* ${RESULTSDIR}/tmp* -fi - -# TODO this now commented out, so to finish inside job folder -#cd - - -# Stop printing commands and their arguments as they are executed. - -echo -echo "*********************************************************************" -echo "* SUCCESS: ${JOB_NAME} PROCESSED WITH ${PIPELINEVERSION} *" -echo "*********************************************************************" -echo -echo Bye! -echo - diff --git a/workflow/config_adapter.py b/workflow/config_adapter.py new file mode 100644 index 0000000..6f57986 --- /dev/null +++ b/workflow/config_adapter.py @@ -0,0 +1,64 @@ +from __future__ import annotations + +from pathlib import Path +from typing import Any + +import yaml + + +# Keys accepted from CWL-style job YAML and merged into Snakemake config. +CWL_COMPAT_KEYS = { + 'dataset', + 'barcode1r1FastqGz', + 'barcode1rmRepBam', + 'barcode1Inputr1FastqGz', + 'barcode1InputrmRepBam', + 'bowtie2_db', + 'bowtie2_prefix', + 'fileListFile1', + 'gencodeGTF', + 'gencodeTableBrowser', + 'repMaskBEDFile', + 'prefixes', + 'se_or_pe', +} + + +def _normalize_value(value: Any, *, base_dir: Path) -> Any: + if isinstance(value, dict): + if 'path' in value and isinstance(value['path'], str): + candidate = Path(value['path']) + if candidate.is_absolute(): + return str(candidate) + return str((base_dir / candidate).resolve()) + return {k: _normalize_value(v, base_dir=base_dir) for k, v in value.items()} + if isinstance(value, list): + return [_normalize_value(v, base_dir=base_dir) for v in value] + return value + + +def load_yaml(path: str | Path) -> dict[str, Any]: + yaml_path = Path(path) + with yaml_path.open('r', encoding='utf-8') as fh: + data = yaml.safe_load(fh) or {} + if not isinstance(data, dict): + raise ValueError(f'Expected mapping at YAML root: {yaml_path}') + return data + + +def map_cwl_job_to_config(job_data: dict[str, Any], *, base_dir: str | Path) -> dict[str, Any]: + based = Path(base_dir) + mapped: dict[str, Any] = {} + for key, value in job_data.items(): + if key in CWL_COMPAT_KEYS: + mapped[key] = _normalize_value(value, base_dir=based) + return mapped + + +def merge_cwl_job_yaml_into_config(config: dict[str, Any], cwl_job_yaml: str | Path) -> dict[str, Any]: + yaml_path = Path(cwl_job_yaml) + job_data = load_yaml(yaml_path) + mapped = map_cwl_job_to_config(job_data, base_dir=yaml_path.parent) + merged = dict(config) + merged.update(mapped) + return merged diff --git a/workflow/rules/se_foundation.smk b/workflow/rules/se_foundation.smk new file mode 100644 index 0000000..a0da631 --- /dev/null +++ b/workflow/rules/se_foundation.smk @@ -0,0 +1,57 @@ +rule all: + input: + 'results/mini/split/verified.ok', + 'results/mini/merge/verified.ok', + + +rule unzip_mini_sam: + input: + config['mini_validation']['source_sam_gz'] + output: + 'results/mini/split/ip.preRmDup.sam.mini.sam' + shell: + 'gzip -cd {input} > {output}' + + +rule split_mini_sam_python: + input: + 'results/mini/split/ip.preRmDup.sam.mini.sam' + output: + directory('results/mini/split/python_tmp') + shell: + r''' + rm -rf {output} + mkdir -p {output} + cd {output} + python3 {workflow.basedir}/bin/python/split_bam_to_subfiles_SEorPE.py ../ip.preRmDup.sam.mini.sam SE + ''' + + +rule verify_split_against_manifest: + input: + split_dir='results/mini/split/python_tmp', + manifest=config['mini_validation']['split_manifest'], + output: + 'results/mini/split/verified.ok' + script: + '../scripts/verify_split_manifest.py' + + +rule merge_parsed_python: + input: + config['mini_validation']['merge_input_1'], + config['mini_validation']['merge_input_2'], + output: + 'results/mini/merge/merged.python.parsed' + shell: + 'python3 {workflow.basedir}/bin/python/merge_multiple_parsed_files.simplified_20191022.py {output} {input[0]} {input[1]}' + + +rule verify_merge_against_expected: + input: + produced='results/mini/merge/merged.python.parsed', + expected=config['mini_validation']['merge_expected'], + output: + 'results/mini/merge/verified.ok' + script: + '../scripts/verify_merge_against_expected.py' diff --git a/workflow/scripts/verify_merge_against_expected.py b/workflow/scripts/verify_merge_against_expected.py new file mode 100644 index 0000000..20ab9e9 --- /dev/null +++ b/workflow/scripts/verify_merge_against_expected.py @@ -0,0 +1,26 @@ +import hashlib +from pathlib import Path + + +def sha256_file(path: Path) -> str: + h = hashlib.sha256() + with path.open('rb') as fh: + while True: + chunk = fh.read(1024 * 1024) + if not chunk: + break + h.update(chunk) + return h.hexdigest() + + +produced = Path(snakemake.input['produced']) +expected = Path(snakemake.input['expected']) +out_ok = Path(snakemake.output[0]) + +assert produced.exists(), f'missing produced file: {produced}' +assert expected.exists(), f'missing expected file: {expected}' +assert produced.stat().st_size == expected.stat().st_size, 'size mismatch' +assert sha256_file(produced) == sha256_file(expected), 'sha mismatch' + +out_ok.parent.mkdir(parents=True, exist_ok=True) +out_ok.write_text('ok\n', encoding='utf-8') diff --git a/workflow/scripts/verify_split_manifest.py b/workflow/scripts/verify_split_manifest.py new file mode 100644 index 0000000..2eaf7bf --- /dev/null +++ b/workflow/scripts/verify_split_manifest.py @@ -0,0 +1,31 @@ +import csv +import hashlib +from pathlib import Path + + +def sha256_file(path: Path) -> str: + h = hashlib.sha256() + with path.open('rb') as fh: + while True: + chunk = fh.read(1024 * 1024) + if not chunk: + break + h.update(chunk) + return h.hexdigest() + + +split_dir = Path(snakemake.input['split_dir']) +manifest = Path(snakemake.input['manifest']) +out_ok = Path(snakemake.output[0]) + +with manifest.open('r', encoding='utf-8') as fh: + rows = list(csv.DictReader(fh, delimiter='\t')) + +for row in rows: + p = split_dir / row['filename'] + assert p.exists(), f'missing split output {p}' + assert str(p.stat().st_size) == row['size_bytes'], f'size mismatch {p}' + assert sha256_file(p) == row['sha256'], f'sha mismatch {p}' + +out_ok.parent.mkdir(parents=True, exist_ok=True) +out_ok.write_text('ok\n', encoding='utf-8')