TE coordinates, TE MSA, and genome-wide dataset are required for the pipeline. Although there are different methods to prepare these files, the format requirement for each input file, and the preparation method and used including parameters and options in each showcase are outlined below:
The requirements for this input is to be in BED-like format (chromosome, start, end, and strand must be present) and optionally provide unique ID for TE entries. It is possible for one ID to have multiple entries in the table as the package supports fragment joining. Critically, the ID from TE coordinate should be the same as TE sequence ID in TE MSA.
For the testing sample, the coordinates for THE1C TEs were obtained from RepeatMasker output. The header rows were edited from two rows into one rows to facilitate importing by pandas. After importing, the THE1C entries were extracted from the table using ma_mapper.sequence_alignment. This module is a wrapper for biopython.SeqIO and MAFFT to streamline the sequence alignment process.
from ma_mapper import sequence_alignment
coord_table=sequence_alignment.extract_coord_from_repeatmasker_table(
subfamily = 'THE1C',
repeatmasker_table = '/path/to/repeatmasker/output/file/',
save_to_file = False
)The resulting coord_table is a pandas.DataFrame object that is in BED format and compatible with the pipeline. To save the table to a file, users can provide a destination path via save_to_file.
The TEA MSA file is required to be in FASTA format with a header specifically generated using their coordinate: >ID::chromosome:start-end(strand) to ensure proper row syncing across matrices in the downstream analyses. For this package, the ma_mapper.sequence_io() was used to extract sequence from the human reference genome (hg38) obtained from UCSC.
from ma_mapper import sequence_alignment
te_seqeunces=sequence_alignment.sequence_io(
coordinate_table=coord_table,
source_fasta='/path/to/hg38/fasta',
save_to_file='/target/path/for/TE/sequences')coordinate_table can accept both a dataframe (as produced in the previous step) and a path to BED file. Although the function returns a SeqIO object, it is recommended to save the result as a FASTA file for use in the next step.
Then, the multiple sequence alignment was performed on TE sequences using ma_mapper.sequence_alignment.mafft_align.
from ma_mapper import sequence_alignment
sequence_alignment.mafft_align(
input_filepath='/path/to/TE/sequences',
nthread=6,
output_filepath='/target/path/for/TE/alignment'
)The output was saved in FASTA format. Users can adjust multithread parameters using nthread, nthreadb, nthreadit and provide additional flags using mafft_arg argument.
To show case the capability of the framework, genome-wide datasets in this thesis were obtained from various publicly available sources as follows:
An example of BED file was obtained from HOMER transcription factor motif prediction. The table was then separated by motif names into smaller table for efficient data retrieval.
An example BAM file was generated from ZNF267 ChIP-exo derived from KZFP overexpression experiments. The FASTQ file, retrieved from Sequence Read Archive (SRA) was mapped to GRCh38_noalt_as reference genome using bowtie2. Subsequently, the mapped reads in SAM format was converted into BAM, then sorted and indexed by samtools. The read mapping operation was performed with default settings as shown below:
# Step 1: Align with Bowtie2
bowtie2 --threads 128 --local -x GRCh38_noalt_as -U /path/to/fastq/ -S /target/path/to/sam/
# Step 2: Convert SAM to BAM
samtools view -@ 128 -bS /target/path/to/sam/ > /target/path/to/bam/
# Step 3: Sort BAM
samtools sort -@ 128 -o /target/path/to/sam/sorted/ /target/path/to/sam/
# Step 4: Index BAM
samtools index /target/path/to/sam/sorted/An example of BigWig file, a phyloP track of Zoonomia multispecies alignment , was obtained from UCSC database and was used without modification.
Example VCF files were obtained from the gnomAD database and stored in the same directory. After downloading, these files were indexed by tabix for efficient data retrieval.
The following snippet is the command for VCF file indexing:
tabix -p vcf /path/to/zipped/vcf/file/\label{mafprep}
Example MAF files of the updated Zoonomia dataset obtained from UCSC database. As MafIO.MafIndex() does not support a compressed MAF file, all files were unzipped before indexing for fast data retrieval. For this framework, .mafindex files are required to have exactly the same name as their .maf counterparts as shown below.
chr1.maf
chr1.mafindex
...
chrY.maf
chrY.mafindexThe MAF files were indexed using the code snippet below:
# build mafindex file
from bio import MafIO
maf_file = '/path/to/maf/'
target_species = 'hg38'
maf_id = f'{target_species}.{chrom}'
mafindex_filedir = '.'.join(str.split(maf_file, sep='.')[:-1]) #remove .maf
mafindex_filepath = f'{mafindex_filedir}.mafindex'
index_maf = MafIO.MafIndex(mafindex_filepath, maf_file, maf_id) Alignment conversion is mandatory for data mapping as the numerical alignment matrix is needed as a reference for gap position. The input for this step is the aligned FASTA file with specific sequence header >ID::chromosome:start-end(strand).
from ma_mapper import mapper
alignment_filepath = '/path/to/alignment/fasta/'
filtered_alignment_matrix, alignment_coordinate = mapper.parse_and_filter(alignment_file=alignment_filepath,col_threshold = 0.10, col_content_threshold = 0.10, row_threshold = 0.10)The outputs of this step are a filtered numerical alignment matrix numpy.ndarray and a metadata table pandas.DataFrame. By default, these outputs are already filtered internally. For pre-filtered output, users can toggle the preprocess_out argument to make the function returns a pre-filtered numerical alignment matrix, a pre-filtered metadata, and filters.
Filtering behavior is customizable with several arguments:
filter: enable/disable filteringrow_threshold: set the row thresholdcol_threshold: set the column thresholdcol_content_threshold: sets the minimum non-gap content for a column
Additionally, to examine regions beyond TE boundaries, users can specify the extension length and the path to the genome sequence for sequence extraction using extension and source_fasta arguments respectively.
The inputs for this step are the numerical alignment matrix, a path to genome-wide data file, and a type of file. Regardless of input type, this function always outputs a numerical data matrix that has the same gapped position as the numerical alignment matrix. The specific handling and parameters for each supported file type is described below.
alignment_filepath = '/path/to/alignment/fasta'
genomewide_data_filepath = '/path/to/bed/file/'
filtered_alignment_matrix, alignment_coordinate = mapper.parse_and_filter(alignment_file=alignment_filepath,col_threshold = 0.10, col_content_threshold = 0.10, row_threshold = 0.10)
bed_data_matrix=mapper.map_and_overlay(alignment_filepath, genomewide_data_filepath,data_format='bed',strand_overlap=True,
col_threshold = 0.10, col_content_threshold = 0.10, row_threshold = 0.10)Some genome-wide dataset BED files may lack strand information (i.e., . instead of +,-). The strand_overlap argument can be used to enable/disable strand-specific processing.
alignment_filepath = '/path/to/alignment/fasta'
genomewide_data_filepath = '/path/to/bam/file/'
filtered_alignment_matrix, alignment_coordinate = mapper.parse_and_filter(alignment_file=alignment_filepath,col_threshold = 0.10, col_content_threshold = 0.10, row_threshold = 0.10)
data_matrix_forward=mapper.map_and_overlay(alignment_filepath, genomewide_data_filepath,data_format='read_forward',
col_threshold = 0.10, col_content_threshold = 0.10, row_threshold = 0.10)
data_matrix_reverse=mapper.map_and_overlay(alignment_filepath, genomewide_data_filepath,data_format='read_reverse',
col_threshold = 0.10, col_content_threshold = 0.10, row_threshold = 0.10)For BAM files, users can choose to separate read_forward and read_reverse or the combined read_sum. Smoothing is also supported via normal_forward and normal_reverse with parameters for for probe length probe_length and smoothing length smoothing_length.
Even though BAM files are recommended to be preprocessed before using this pipeline, basic filtering using pysam internal flags is supported. Users can enable arguments such as:
drop_duplicate: remove duplicate readsdrop_unmapped: remove unmapped readsdrop_2: remove secondary/supplementary readsdrop_multipmap: remove multimapped reads
All filters are disabled by default.
alignment_filepath = '/path/to/alignment/fasta'
genomewide_data_filepath = '/path/to/bigwig/file/'
filtered_alignment_matrix, alignment_coordinate = mapper.parse_and_filter(alignment_file=alignment_filepath,col_threshold = 0.10, col_content_threshold = 0.10, row_threshold = 0.10)
bigwig_data_matrix=mapper.map_and_overlay(alignment_filepath, genomewide_data_filepath,data_format='bigwig',
col_threshold = 0.10, col_content_threshold = 0.10, row_threshold = 0.10)There is no additional parameter for BigWig.
alignment_filepath = '/path/to/alignment/fasta'
genomewide_data_path = '/path/to/vcf/directory/'
vcf_data_matrix=mapper.map_and_overlay(alignment_filepath, genomewide_data_path,data_format='vcf',query_key='AF',vcf_format='gnomad',
col_threshold = 0.10, col_content_threshold = 0.10, row_threshold = 0.10)Here, vcf_format argument has two options, one for gnomAD dataset (fully supported) and another for direct file extraction. For gnomAD extraction, users have to provide directory path to gnomAD VCF files. Other than alternate allele frequency value, query_key argument also allows the users to access other fields.
alignment_filepath = '/path/to/alignment/fasta'
genomewide_data_path = '/path/to/maf/directory/'
maf_data_matrix=mapper.map_and_overlay(alignment_filepath, genomewide_data_path,data_format='maf',separated_maf=True, count_arg='common_freq', target_species='hg38',
col_threshold = 0.10, col_content_threshold = 0.10, row_threshold = 0.10)Similarly, separated_maf toggle allows either per-chromosome MAF files or a single file extraction. For separated MAF extraction, users must provide directory path. The target_species argument sets the reference genome, for example, hg38 for the human genome in the updated Zoonomia dataset.
The count_arg parameter control the format at nucleotide level. Avialable options are:
ref_freq: reference base frequencycoverage: number of species aligned (not to be confused with sequencing coverage)common_freq: most common base frequencycommon_freq_nogap: same as above, excluding gaps from total countcommon_raw: raw count of the most common basebase_count: all base countsbase_freq: all base frequenciesbase_freq_nogap: all base frequencies excluding gapsraw: returns apandas.DataFramewith alignment content at each positionraw_genomesame asrawbut only with genome names only, excluding chromosome labels.
The plots module was used to visualize the resulting data matrix from the previous step. The code below is the intended simplest representation of the matrix.
from ma_mapper import plots
from ma_mapper import mapper
plots.plot(
data = [bigwig_data_matrix],
alignment=filtered_alignment_matrix,
show_alignment=False,
heatmap_color=['RdBu'],
heatmap_mode='overlay',
heatmap_title=['Example of plot function'],
heatmap_title_fs = 10,
heatmap_xlabel = 'position (bp)',
heatmap_xlabel_fs = 10,
heatmap_ylabel = 'sequences',
heatmap_ylabel_fs = 10,
vlim = [[-1,1]],
colorbar=True,
colorbar_steps=[0.1]
)
Slicing/Cropping for specific parts of the plot is possible using plot coordinates, as shown below:
from ma_mapper import plots
from ma_mapper import mapper
plots.plot(
data = [bigwig_data_matrix],
alignment=filtered_alignment_matrix,
show_alignment=False,
heatmap_color=['RdBu'],
heatmap_mode='overlay',
heatmap_title=['Example of plot function'],
heatmap_title_fs = 10,
heatmap_xlabel = 'position (bp)',
heatmap_xlabel_fs = 10,
heatmap_ylabel = 'sequences',
heatmap_ylabel_fs = 10,
vlim = [[-1,1]],
colorbar=True,
colorbar_steps=[0.1],
xlim=[200,250],
)
Additionally, the metadata table which contains TE coordinates, can be expanded further with additional information from either external sources or other data matrices transfromed into row-wise data arrays and subsequently converted into plot annotation
See the tutorials here for more complex usages such as multiple matrix overlays or row annotation.