Motif prediction of abemaciclib in a breast cancer cell line using ChIP-Seq data analysis

Abstract

Introduction: Chromatin immunoprecipitation sequencing (ChIP-Seq) is a DNA sequencing technique for the identification of binding sites in genomic sequences. ChIP-Seq experiments are a combination of immunoprecipitation and sequencing techniques that are used for the identification of chromatin regions that bind molecules such as transcription factors (TFs), histones, and drugs. In this study, computational analysis of ChIP-Seq data was performed to predict the binding sites in breast cancer cells and their association with several molecules, such as TFs and drugs. A complete and comprehensive computational study has been performed to predict the binding sites of abemaciclib. Functional enrichment of selected motifs was performed to identify important motifs that function in breast cancer and show binding with the drug abemaciclib.

Methods and Materials: The ChIP-Seq analysis protocol was performed using the Galaxy server (https://usegalaxy.org/). The abemaciclib binding motif was identified using MEME tools. For this research, ChIP-Seq data from a breast cancer cell line was retrieved from the GEO database, accession number GSM4763932. This dataset includes ChIP-Seq of MCF-7 cells exposed to abemaciclib. The ENA browser was used to retrieve the data. Statistical analysis was performed using the default parameters of Fast-QC, Multi-QC, Map with BWA, Filter SAM and BAM, MACS2, and ChIP Seeker tools on the Galaxy server.

Results: Computational analysis identifies the abemaciclib consensus binding sequence as TGGCTCACGCCTGTAATCCCAGCACTTT, and this motif occurs 2980 times in the Homo sapiens reference genome hg19.

Conclusions: This study identifies the binding sites and affinity of abemaciclib in a breast cancer cell line.

Introduction

Cancer of the breast is uncontrolled production of breast cells1. Breast cancer is the second most common cancer diagnosed in women in the United States after skin cancer. Breast cancer may occur in both men and women, but it is much more prevalent in women. It is the most prevalent cause of death from cancer among women worldwide2. In more developed countries, incidence rates are high, while rates in less developed countries are low but rising rapidly.

Genomics studies have revealed genes associated with breast cancer. The identification and functional understanding of these genes are very important for the molecular diagnostics and drug development industries. Breast cancer genes such as BRCA1, BRCA2, CHEK2, ATM and PALB2 have been widely studied since these genes play a major role in breast cancer development3. Cyclin-dependent kinases 4 and 6 (CDK4/6) are pharmacological inhibitors and have a significant effect on the practice of oncology4. They are routinely recommended for estrogen receptor positive breast cancer therapy, and many studies are underway to assess their activity against other forms of cancer5. CDK4/6 mediates the cell cycle transition from G1 to S phase, and CDK4/6 inhibitors cause arrest of the G1 phase in tumor cells6. Other phenotypes in cancer cells could also be activated, including enhanced immunogenicity, apoptotic evasion, histologic tumor differentiation, and increased dependence on receptor tyrosine kinase signaling7. Abemaciclib is a CDK4- and CDK6-selective ATP-competitive, reversible kinase inhibitor that has demonstrated antitumor activity in clinical trials as a single agent in hormone receptor-positive (HR+) metastatic breast cancer8.

With advancements in next-generation sequencing technologies, genome-wide screening of binding sites can be performed by sequencing chromatin associated with proteins. Chromatin immunoprecipitation (ChIP-Seq) is used to identify the binding sites of proteins of interest at the genomic scale. ChIP-Seq analysis is used to identify the binding sequence for any drug or transcription factor9. Transcription factors are proteins that bind upstream of genes to enhance or initiate transcription10. Many genes are responsible for causing breast cancer, as described above11. As the main genes responsible for tumor formation in breast cancer are BRCA1 and BRCA2, ChIP-Seq can be used to identify other transcription factors that initiate the activation of BRCA1 or BRCA2 and facilitate the formation of tumors, causing cancer12.

ChIP-Seq is an effective technique for defining genome-wide DNA-binding sites for transcription factors and other proteins by combining chromatin immunoprecipitation (ChIP) with sequencing13. The DNA-bound protein is immunoprecipitated using a specific antibody14. Bound DNA is co-precipitated, purified, and sequenced. ChIP-Seq can be used to identify the binding sites of proteins associated with DNA and can be used to map a given protein's genome-wide binding sites15. Usually, ChIP-Seq begins with DNA–protein complex crosslinking. The samples are fragmented and processed to trim unbound oligonucleotides with an exonuclease. To immunoprecipitate a DNA–protein complex, protein-specific antibodies are used16. The DNA is extracted and sequenced, thereby identifying the sequences of the protein-binding sites with high resolution17. Insights into gene regulation events that play a role in different diseases and biological pathways, such as development and cancer progression, have been uncovered by the application of next-generation sequencing (NGS) to ChIP. ChIP-Seq allows the interactions between proteins and nucleic acids to be thoroughly studied on a genome-wide scale18.

Michigan Cancer Foundation-7 (MCF-7) is a breast cancer cell line that was first isolated from a 69-year-old Caucasian women in 197019. It is a well-known cell line used worldwide. It has features of mammary epithelial differentiation and is positive for many epithelial markers, such as β-catenin, cytokeratin 18 (CK18) and E-cadherin20. It is also negative signal for mesenchymal markers such as smooth muscle actin and vimentin21. MCF‑7 cells maintain the expression of other specific molecular markers of epithelium and do not express CD-44 ligand.

Abemaciclib is an antitumor agent and a dual cyclin-dependent kinase 4 (CDK4) and 6 (CDK6) inhibitor. CDK4/6 are involved in the cell cycle and in the event of unregulated activity, promote cancer development22. On September 28, 2017, the FDA approved abemaciclib for the treatment of HR-positive and HER2-negative advanced or metastatic breast cancer that has progressed following failed endocrine therapy, under the brand name Verzenio23. It is an oral inhibitor of cyclin-dependent kinase (CDK) that targets the cell cycle pathways of CDK4 (cyclin D1) and CDK6 (cyclin D3) with potential antineoplastic activity24. Abemaciclib inhibits CDK4 and 6 directly, thus inhibiting the phosphorylation of retinoblastoma (Rb) protein in early G125. Inhibition of Rb phosphorylation prevents the transition of the CDK-mediated G1-S process, thereby arresting the cell cycle in G1 phase, suppressing DNA synthesis and inhibiting the growth of cancer cells26. The interaction of abemaciclib with different genes such as AKT1, BARD1, and CDK4, and the functional consequences thereof, is shown in Table 1. Identification of binding sites for drugs is important to understand their function. Any mutation in the genome of the cell that disrupts a binding sites will result in loss of binding of the drug. In order for a drug to be effective in a mutated cell, such as a cancer cell, it is essential for the drug to bind its target sequence. Keeping this in mind, breast cancer data were used in the current study to identify the core binding sites in the genome of a breast cancer cell line. Functional analysis was also performed to identify functional sequences and core binding sites.

Table 1.

Genes accosiated with Abemaciclib Drug

S.No	PubChem Gene	Interaction with Abemaciclib drug	Evidence PMID
1	ABCB1 and ABCB2	Abemaciclib results in decreasing the activity of ABCB1 and ABCB2 protein 27	27816545
2	AKT1	Abemaciclib results in decreased phosphorylation of and results in decreased activity of AKT1 protein 28	26909611
3	BARD1	Abemaciclib results in decreased expression of BARD1 mRNA 29	28620137
4	CDK4	Abemaciclib results in decreased activity of [CDK4 protein binds to CCND1 protein] 30	24919854
5	RB1	Abemaciclib results in decreased phosphorylation of RB1 protein 31	26909611
6	SPARC	Abemaciclib results in increased expression of SPARC mRNA 32	28620137
7	WNT4	Abemaciclib results in decreased expression of WNT4 mRNA 33	28620137

Methods

Data Retrieval

The SRA (Short Read Archive) database was used for ChIP-Seq data retrieval. GSM4763932 ID was seen relevant because the study was on the MCF-7-cell line treated with abemaciclib. The SRA ID SRA1120692 was used. Two samples with accession numbers SRR12576544 and SRR12576545 were selected for motif enrichment and analysis. The ENA browser was used to retrieve the data (https://www.ebi.ac.uk/ena/browser/text-search?query=SRA1120692). The SRA IDs of MCF-7 DMSO and MCF-7 abemaciclib samples are shown in Table 2. FASTQ files were uploaded to the Galaxy server https://usegalaxy.org/ for analysis and annotation.

The Galaxy (https://usegalaxy.org/) web platform was used for the analysis and annotations of ChIP-Seq data. The methodology is shown in Figure 1. The tools used in Galaxy were Fast-QC, Multi-QC, Map with BWA, Filter SAM and BAM, MACS2, and ChIP Seeker. MEME Multiple EM for motif elicitation was used for motif identification at http://meme-suite.org/tools/meme. The steps followed for the identification of motifs are described below.

Step I: Fast-QC34

The Fast-QC tool was used to perform a quality check on raw sequence data. This offers a modular collection of tests that can be used to provide an indication of whether data have any issues that should be resolved before further processing. SRR12576544 and SRR12576545 were used as inputs for the Fast-QC analysis. The outputs were in HTML and .txt format. Fast‑QC tools produce statistical values that can be interpreted to identify filtering or preprocessing steps.

Step II: Multi-QC35

Multi-QC aggregates data from different samples into a single report. It compiles all results into a single HTML report that can be used for comparative study. The inputs for the Multi‑QC tool were the Fast-QC reports, and the compiled results were in the form of plots.

Step III: Map with BWA36

This tool is used to align short sequences to the large sequence database. SRR12576544 and SRR12576545 in fasta.qz file format were used as inputs for the map with BWA tool. The output was in BAM format, which can be visualized in the UCSC, bam.iobib.io and IGV genome browsers.

Step IV: Filter SAM or BAM37

This tool uses the samtools view commands in SAM Tools to filter SAM and BAM files on the mapping quality. This is a unique way of mapping by removing the non-required sequence. The input file is the result of the map with the bam tool, and the output is in BAM format.

Step V: Multi Bam Summary38

For two or more BAM files, multiBamSummary calculates the read coverage for genomic regions. The input was the result from the Filter SAM or BAM tool, and output was in the deeptools_coverage_matrix format. This result can be used for comparison of genomic alignments in different samples.

Step VI: plot Correlation39

This tool was used to analyze and visualize sample correlations based on multiBamSummary or multiBigwigSummary. The input was the result of multiBamSummary, and the output was in png format. This plot was used for the comparison of variations in sample files.

Step VII: plot fingerprint40

This tool samples indexed BAM files and graphs a cumulative read coverage profile for each file. The input was the same as for plot correlation, and the output was also in png format. This graph gives detailed information about the location of binding sites in genomic regions.

Step VIII: bamCoverage41

This method takes an alignment of reads or fragments as input (BAM file) and produces as output a coverage track (bigWig or bedGraph).

Step IX: MergeSamFiles42

Similar to the "merge" feature of Samtools. This tool is used to combine SAM or BAM files from separate runs and read groups into a single file.

Step X: Split 43

This is used to separate a combined file from a single BAM file. Aligned reads corresponding to the original four datasets are included in every subsequent BAM file.

Step XI: MACS244

ChIP-Seq model-based analysis (MACS) is a widely used method for identification of transcription factor binding sites. To determine the importance of enriched ChIP regions, the MACS algorithm captures the effect of genome complexity. It is also suitable for wider areas, although it was developed for the detection of transcription factor binding sites. The input was the result of Split, and the output was in the form of bed format for submit peaks and narrow peaks.

Step XII: ChIPSeeker45

ChIPseeker is a Bioconductor package in Galaxy for annotating ChIP-Seq data. The annotatePeak feature performs peak annotation. In addition to the distance from the peak to the TSS of its closest gene, the location and strand information of the nearest genes were recorded. The input file is obtained from the gencode tool, and the output is in the form of plots.

Step XIII: MEME-ChIP 46

MEME-ChIP performs motif analysis on large sets of sequences, such as those identified by ChIP-Seq experiments. The input was obtained from the Extract Genomic DNA tool, and the results were in HTML format.

Step XIV: plot Heatmap 47

This tool creates a heatmap for scores associated with genomic regions. This plot is used to visualize the genomic locations of predicted binding sites.

Download figure
Enlarge figure

Figure 1 . Workflow used for ChiP-Seq data analysis and motif identification . This workflow shows the data pre-processing steps, tools used for genomic alignment, prediction of motifs, and functional enrichment of motifs.

Download figure
Enlarge figure

Figure 2 . Number of duplicate and unique reads in ChiP-Seq files of MCF-7 DMSO (SRR12576544) and MCF-7 Abemaciclib (SRR12576545) . This comparison shows the number of duplicated reads in both the samples, that is 15 M in (SRR12576544) and 20 M in (SRR12576545).

Download figure
Enlarge figure

Figure 3 . Genomic alignment of SRR12576544 to reference genome hg19 . This data details the number of read alignments in each chromosome and for the visualization of genomic alignment. The result shows 92.6% mapped reads across the genome.

Download figure
Enlarge figure

Figure 4 . Genomic alignment SRR12576545 to reference genome hg19 . The result shows 96 % mapped reads across the genome, that is much higher than compared to SRR12576544.

Download figure
Enlarge figure

Figure 5 . Heatmap Correlation between aligned reads of both the files MCF-7 DMSO (SRR12576544) and MCF-7 Abemaciclib (SRR12576545) . This data compares the percentage of aligned reads in approximately same sample file.

Download figure
Enlarge figure

Figure 6 . Fingerprint plot of both the samples produced using plot fingerprint tool . It samples indexed BAM files and graphs a cumulative read coverage profile for each sample. It shows (SRR12576545, red line) has more read coverage.

Download figure
Enlarge figure

Figure 7 . MACS2 peak result shows the percentage of DNA-Protein interaction in both the strands forward and reverse strands . Enrichment of binding sites in both the strands is compared to verify the predicted binding sites.

Download figure
Enlarge figure

Figure 8 . Pie graph of identified peaks in MCF-7 Abemaciclib (SRR12576545) . Across different genomic functional regions with percentages identified using ChIPSeeker. It shows binding sites are present in functional regions like promoter 11.82%, Exon 0.07%, and others are in the intergenic region 27.2%.

Download figure
Enlarge figure

Figure 9 . MEME motifs found in 2980 sequences corresponding to common peak regions in MCF-7 Abemaciclib (SRR12576545) . This motif was predicted using the MEME ChIP tool, and it shows that the Abemaciclib drug has a binding affinity toward this predicted motif.

Download figure
Enlarge figure

Figure 10 . Identification of motif “TGGCTCACGCCTGTAATCCCAGCACTTT”, which was identified in the 2980 positions in the Homo sapiens reference genome version hg19 . This details about presence of predicted motif in different chromosomes along with their start and end positions.

Download figure
Enlarge figure

Figure 11 . Heatmap plot shows the percentage of predicted motifs across the genome . Blue region in the plot depicts the region of unique motifs predicted in the sample files. Heatmap results confirms that Abemaciclib has binding regions in this breast cancer cell line and predicted motifs have function in cells.

Table 2.

SRR ids with their GSM number used for this study

S.No	Name	SRR number	GSM number
1	MCF-7 DMSO	SRR12576544	GSM4763931
2	MCF-7 Abemaciclib	SRR12576545	GSM4763932

Table 3.

Multi-QC result showing Duplicates, GC content and Maximum length sequence

Sample Name	% Dups	% GC	M Seqs
SRR12576544_fastq_gz	31.0%	49%	23.6
SRR12576545_fastq_gz	8.0%	46%	21.8

Results

Following the protocol described in Figure 1, Fast-QC was performed for SRR12576544 and SRR12576545 samples in Galaxy, followed by Multi-QC. The input file used for this workflow was in fastq format, and the output file generated from Multi-QC was in HTML format. The results of Multi-QC are shown in Table 3. The comparison between sample files in terms of the number of duplicate reads, GC content and maximum length (M Seq) shows that both files can be used for further study. Figure 2 shows the number of duplicate and unique reads in both samples.

The BAM file was generated from the fastq.gz file by using Map with the BWA tool of Galaxy. The alignment between the sample file and the reference genome, hg19, was performed with the tool Map with BWA to identify the similarity between read sequences and genomic sequences. The results were in bam.iobio, IGV and UCSC formats, which were visualized using https://bam.iobio.io/home, Integrative Genomics Viewer (IGV) (https://software.broadinstitute.org/software/igv/) and the UCSC Genome Browser (https://genome.ucsc.edu/), respectively. The data used for this study satisfied all the parameters of quality control and were of very high quality, so the data were not trimmed or postprocessed in any way before mapping. Figure 3 and Figure 4 show the result of alignment of the read sequences with the hg19 reference genome. Figure 3 shows the alignment between the MCF-7 vehicle control reads and the reference human genome (SRR12576544 to hg19). Figure 4 shows the alignment between the MCF-7 abemaciclib sample and the reference genome (SRR12576545 to hg19). Both samples had more than 95% alignment with the reference genome.

All non-uniquely mapped reads were identified. This was accomplished by simply filtering out all reads with a mapping quality of less than 20 using NGS:SAMtools (SAM or BAM filter). Then, the mapping quality of the mapped reads was analyzed. The multiBamSummary tool was used to calculate the average read coverage for a list of two Bam files. Then, the plotCorrelation tool was used to create a heatmap of correlation scores between the samples, as shown in Figure 5, demonstrating that both files have approximately the same number of aligned regions.

plotFingerprint was used to identify the signal strength of the ChIP sample. Figure 6 shows the comparison of MCF-7 DMSO (SRR12576544, red) and MCF-7 abemaciclib (SRR12576545, green), which shows that there might be some variation in the aligned regions that can be explored in further studies.

The input for the MACS2 peak calling tool was obtained by merging the BAM file with the tool MergeSamFiles, and then the split tool was used to separate the merged file into the individual BAM files. Aligned reads corresponding to the original two datasets were included in each subsequent BAM file. Peak calling was performed using the MACS2 package, the results of which are shown in Figure 7. This tool finds the best parameters for the MACS2 peak calling.

The ChIPseeker tool was used to annotate the ChIP-Seq data. It annotates ChIP peaks and provides functionality to visualize the coverage of ChIP peaks over chromosomes and peak profiles that bind to TSS regions. ChIP peak profile comparison and annotation are also provided. In addition, it supports the evaluation of substantial overlap between ChIP-Seq datasets. Currently, ChIPseeker includes 17,000 GEO database bed file details. The input was given after performing the gene pool, which helped identify promoter, exon, intron, UTR, and downstream regions.

The output was in the form of plots and Excel files, which give all introns and exons with chromosome numbers. The result of the peak annotation is shown in Figure 8, which shows the percentage of identified peaks in promoter, exon, intron, and downstream regions. This is important to identify functional regions that show alignment with sample files. The reference genome result shows that there are alignment regions in almost all functional regions of the genome, such as promoters, genes, and exons.

De novo motif discovery, motif enrichment analysis, motif position analysis and motif clustering were performed using MEME-ChIP, providing a detailed image of the motifs that are enriched in the input sequences. Two complementary types of de novo motif discovery are performed by MEME-ChIP: weight matrix-based discovery for high accuracy and word-based discovery for high sensitivity. The study of motif enrichment using human, mouse, worm, fly and other model organism motifs offers even greater sensitivity.

The Extract Genomic DNA tool was used to generate genomic sequences corresponding to the ChIP peaks. Figure 9 shows the motif that was predicted using MEME ChIP, and Figure 10 shows the motif TGGCTCACGCCTGTAATCCCAGCACTTT, which was identified in 2980 positions in the Homo sapiens reference genome hg19.

A heatmap was generated to visualize the position of the motif in the whole genome (Figure 11). The blue region in the plot depicts the region of all motifs in the whole genome, showing the regions across the genome that represent motifs and binding sites. Heatmap results verify that abemaciclib has binding regions in this breast cancer cell line, and motifs identified are true motifs that have function in cells.

Discussion

ChIP-Seq is a next-generation sequencing technology that is used to identify binding sites genome-wide. This technique has been explored in the drug development industry to characterize the effects of drugs and identify conserved regions of their binding sites. Since drug binding site identification is important to study the functionality and effectiveness of drugs, ChIP-Seq is used to study novel drugs and their effectiveness. It can also be used to screen ligand libraries for binding efficiency in cancerous cells48. Studies of various diseases and cancers have identified variations in binding sites and differential binding sites for drugs49. ChIP-Seq technology has been explored in different dimensions to study breast cancer cells and to identify novel drugs that can stop breast cancer50. Studies have been carried out mainly to identify drugs against breast cancer51. In this study, computational analysis of ChIP‑Seq data of breast cancer cells was performed to identify binding sites of the drug abemaciclib. The binding motif for abemaciclib was identified as TGGCTCACGCCTGTAATCCCAGCACTTT. Mapping of the identified motif to the human reference genome shows that the predicted motif occurs in 2980 positions in the hg19 reference genome. The results identify the binding sites for abemaciclib in the MCF-7 genome. This study can be used to better understand the function of abemaciclib against breast cancer and provide insight into its conserved and functional region in the human genome.

Conclusions

Breast cancer is the most common type of cancer in women. Breast cancers can become metastatic when breast cancer cells spread to different organs of the body through blood vessels or nodes. Currently, abemaciclib is used to treat breast cancer, and studies have been performed to identify the role of abemaciclib in breast cancer. In this research, ChIP-Seq data were retrieved from the SRA database, and the interaction between abemaciclib and chromatin was studied. These data were selected to identify the motif that abemaciclib binds to, which can be explored in different dimensions to understand the function of the motif and the drug.It can be used in toxicology studies, or in studies of the effect of abemaciclib in treating breast cancer. Further verification and wet lab study are required to study the function of the identified motif and its interaction with abemaciclib.

Abbreviations

ChIP-Seq: Chromatin immunoprecipitation sequencing, TFs: transcription factors, MCF: Michigan Cancer Foundation-7, SRA: Short Read Archive, QC: Quality Control, SAM: Sequence Alignment Map, BAM: Binary Alignment Map, MACS2: Model-based Analysis of ChIP-Seq, MEME: Multiple Expectation maximizations for Motif Elicitation

Acknowledgments

We would like to acknowledge Amity Institute of biotechnology, Amity University Uttar Pradesh, Lucknow campus for providing us facilities to conducting this study. This research project is not funded by any specific grant from funding agencies in the public, commercial, or non-profit sectors.

Author’s contributions

Author Utkarsh Deep has done the work on methodology and interpretation of results and Ruchi Yadav contributed as framing methodology and analysis of the results. All authors read and approved the final manuscript.

Funding

None.

Availability of data and materials

Data and materials used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.