cv_coverage

This code compares the coverages of a k-mer group around crosslink events (landmarks) across multiple CLIP datasets. The comparison is conducted within user-specified genomic regions.

Author: aram.amalietti@gmail.com

Dependencies (these are the versions the script was developed with, pandas >= 1 introduced breaking changes, please use these versions):

python=3.7  
pandas=0.24.2  
numpy=1.19.2  
pybedtools=0.8.1  
matplotlib=3.3.2
seaborn=0.11.0

Usage:

python3 <path_to_script> <xl_in> <motifs> <regions> <kmer_len> <fasta> <fai> <regions_file> <smothing> <percentile> <window> <use_scores> <n_cores> <chunk_size> <cap>

• xl_in list of BED files containing landmarks around which to display the motifs, given as path_to_file1,path_to_file2,path_to_file3 .... It can be a single file or a list of files, the plots will be adjusted accordingly;
• motifs group of motifs to be displayed around the input landmarks, for example AAAA,CCCC,GGGG ...;
• regions group of regions to be considered for analysis. Available regions are genome (intron, CDS, UTR3, UTR5, ncRNA, intergenic), whole_gene (intron, CDS, UTR3, UTR5), intergenic, intron, ncRNA, other_exon (UTR5, CDS), UTR3 and UTR5. Single or multiple regions can be analysed and the plots will adjust accordingly. For example, intron,ncRNA;
• kmer_len the length of analysed motifs (in bases);
• fasta is the path to the genome in fasta format;
• fai is the path to the genome index file;
• regions_file is a segmentation file in GTF format, usually obtained with iCount segment but can also come from other sources;
• smoothing is the size of the smoothing window, usually 12;
• percentile defines the percentile for thresholding genomic landmarks by score, if percentile is set to None no thresholding is used and all landmarks will be used for the analisys;
• window flanking distance around the landmarks that is analysed and displayed;
• use_scores If True, each landmark is weighted by its score, else all landmarks will have equal weigths of 1;
• n_cores is the number of threads used in the process (4 is the usual value);
• chunk_size is the number of rows per thread (10000 is the usual value);
• cap is the max value for any landmarks score. It is only used if use_scores is set to True (20 is the recommended value);

Common issues

The script needs writing permission in the staging directory to make results directory and environment variable TMPDIR for temporary files. If you get KeyError: 'TMPDIR' a solution would be to type export TMPDIR=<path_to_folder> in terminal where you want to run the script.

Outputs

• A pdf file with graphs showing % k-mer group coverage around landmarks for each analyzed CLIP dataset;
• A tsv file with % coverage values (raw values used for plotting);
• A text file that saves the number of analyzed landmarks in each sample;
• A tsv file with saved run parameters.