A method for predicting new antigens based on peptide structural indicators

Abstract

The present invention discloses a method for predicting new antigens based on polypeptide structural indicators, which belongs to the field of bioinformatics technology. The method comprises the following steps: S01, obtaining sample sequencing data; S02, identifying sample somatic cell variation; S03, filtering sample variation; S04, translating variation sequences into proteins; S05, segmenting protein fragments, and filtering polypeptides in the normal human proteome; S06, predicting protein fragment structures; S07, calculating corresponding polypeptide structural indicators; S08, screening new antigens. The invention is based on indicators of polypeptide structure, which can characterize the properties of new antigens from a microscopic level, and is an important supplement to the method of screening new antigens based on molecular biological characteristics.

Description

Technical Field

[0001] The present invention belongs to the technical field of tumor immunotherapy bioinformatics, and specifically relates to a method for predicting new antigens based on polypeptide structural indicators. Background Art

[0002] Malignant tumors are one of the most serious diseases that harm humans. Treatments for malignant tumors have made great progress in the past few decades. At present, conventional treatments for malignant tumors include surgery, radiotherapy, chemotherapy, and targeted therapy. However, the above treatments have certain limitations and are easily affected by toxic side effects and tumor recurrence. In recent years, immunotherapy based on activating the immune system to inhibit and kill tumor cells has become a new hot spot in the field of malignant tumors. Among them, immunotherapy based on new antigens has the characteristics of a wide range of effects and few toxic side effects.

[0003] Neoantigen vaccines have significant tumor immunotherapy effects, can be applied to multiple cancer types, and have relatively few toxic side effects. They have become an important member of the immunotherapy family. The core of this treatment method is the selection of immunogenic antigenic peptides. Currently, the prediction of neoantigens includes analyzing whole exome sequencing and transcriptome sequencing data of tumors and normal tissues, identifying DNA mutations in protein coding regions and human leukocyte antigen subtypes, and using bioinformatics methods to obtain mutant peptides translated from mutant DNA, and finally predicting whether the mutant peptides can be presented to the cell surface by human leukocyte antigens.

[0004] The current basis for screening potential new antigens is mainly based on molecular biological characteristics such as mutation abundance, expression level and affinity between the peptide and the corresponding HLA. However, more and more experimental results show that there are factors independent of these characteristics that also affect the immunogenicity of the peptide. For peptides bound to MHC proteins, the influence of characteristics such as charge and hydrophobic surface is determined by the conformation of the peptide in the binding groove and the size and position of various amino acid side chains. Therefore, the prediction of immunogenic peptides can be enhanced by considering the structural characteristics of the peptide / MHC complex. Therefore, in order to improve the screening efficiency of new antigens and the accuracy of predicting the immunogenicity of peptides, it is of great significance to incorporate peptide structural indicators into the peptide immunogenicity screening method. Summary of the invention

[0005] In response to the above-mentioned problems in screening new tumor antigens, the present invention fully considers the influence of polypeptide structure on polypeptide immunogenicity and develops a bioinformatics method for screening new antigens based on polypeptide structural indicators.

[0006] The present invention discloses a method for predicting new antigens based on polypeptide structural indicators, comprising the following steps:

[0007] S01, obtaining sample sequencing data;

[0008] S02, identification of sample somatic variants;

[0009] S03, sample variation filtering;

[0010] S04, translation of variant sequence into protein;

[0011] S05, protein fragment segmentation and filtering of peptides in the normal human proteome;

[0012] S06, predicting protein fragment structure;

[0013] S07, calculating the structural index of the corresponding polypeptide;

[0014] S08, screening for new antigens.

[0015] Preferably, the sample is a fresh tumor tissue sample; as an alternative, a peripheral blood sample can be selected.

[0016] In some embodiments of the present invention, in S01, DNA sequencing is performed on a tumor sample and a normal control sample to obtain sample sequencing data, wherein the sample sequencing data includes a plurality of overlapping or partially overlapping short read sequences.

[0017] In some preferred embodiments of the present invention, in S01, no less than 500X deep sequencing data of tumor samples and 200X deep sequencing data of normal controls are obtained.

[0018] In some preferred embodiments of the present invention, in S01, whole exome sequencing or panel capture sequencing is used.

[0019] In some preferred embodiments of the present invention, S02 includes the following steps:

[0020] S21, aligning the sequencing data of the tumor sample and the normal control sample to the reference genome, and locating the position of the short read sequence on the reference genome;

[0021] S22, removing PCR duplication from the comparison result file;

[0022] S23, identification of sample cell variation;

[0023] Preferably, the method further includes filtering the sample sequencing data, removing low-quality reads and reads containing adapter or sequencing primer sequences, and performing base quality correction on the comparison results.

[0024] In some embodiments of the present invention, S03 includes the following steps:

[0025] S31, filtering out variants with poor credibility;

[0026] S32, filtering variants located in intergenic regions and introns as well as synonymous variants based on the annotation results;

[0027] Preferably, the variants with poor credibility include variants with fewer supporting reads and also appearing in normal control samples.

[0028] In some embodiments of the present invention, in S04, among the variant results that can produce protein sequence changes after filtering, a mutant transcript is constructed based on the genome mutation information and annotation information and translated according to the translation rules to obtain a mutant protein sequence.

[0029] In some embodiments of the present invention, in S05, the obtained protein sequence is divided into Kmers with a length of 9 to 12, and polypeptides in the normal human proteome are filtered.

[0030] In some embodiments of the present invention, in S06, the segmented protein sequence is input using the software pep-fold, and according to the corresponding algorithm, the tertiary structure diagram of the polypeptide and the protein three-dimensional structure data file can be obtained.

[0031] In some embodiments of the present invention, in S07, the result file of the structural index is obtained using Gromacs software, and the immunogenic polypeptides are screened out through analysis.

[0032] In some embodiments of the present invention, S08 includes the following steps:

[0033] S81, preliminary screening of peptides using RMSF;

[0034] S82, further screening was performed using SASA;

[0035] Preferably, the difference in screening SASA is between 0.2 and 1 nm 2 between the peptides.

[0036] Compared with the prior art, the solution of the present invention has the following advantages:

[0037] The present invention is based on indicators of polypeptide structure and can characterize the properties of new antigens at a microscopic level, thus being an important supplement to the method of screening new antigens based on molecular biological characteristics. BRIEF DESCRIPTION OF THE DRAWINGS

[0038] Figure 1 A flowchart of a method for predicting new antigens based on polypeptide structure indicators according to one embodiment of the present invention. DETAILED DESCRIPTION

[0039] The following is an explanation of the embodiments of the present invention by specific examples, and those skilled in the art can easily understand other advantages and effects of the present invention from the contents disclosed in this specification. The present invention can also be implemented or applied through other different specific embodiments, and the details in this specification can also be modified or changed in various ways based on different viewpoints and applications without departing from the spirit of the present invention. Unless otherwise specified, the embodiments and comparative examples are parallel tests with the same components, component contents, preparation steps, and preparation parameters.

[0040] Example 1

[0041] like Figure 1 As shown, a method for predicting new antigens based on peptide structure indicators:

[0042] S01, obtain sample sequencing data

[0043] Specifically, DNA sequencing is performed on tumor samples and normal control samples to obtain sample sequencing data. The sample sequencing data includes multiple overlapping or partially overlapping short read sequences. The degree of overlap is related to the depth of sequencing. No less than 500X deep sequencing data of tumor samples and 200X deep sequencing data of normal controls should be obtained.

[0044] Preferably, the present invention can use whole exome sequencing or panel capture sequencing.

[0045] S02, Identify sample somatic variants

[0046] Specifically, the sequencing data of tumor samples and normal control samples are first aligned to the reference genome, respectively, to locate the position of the short read sequence on the reference genome. Specifically, BWA software can be used to perform sequence alignment on the sample filtered data.

[0047] Secondly, the comparison result file is processed to remove PCR duplicates. In the process of preparing the library, due to some deviations in the PCR amplification process, some sequences will be over-amplified. These over-amplified reads cannot be used as evidence for variant detection. Therefore, these duplicates formed by PCR amplification should be removed as much as possible. Specifically, the software Picard can be used to compare the file to remove PCR duplicates.

[0048] Thirdly, identify the cell variation of the sample. Somatic cell variation refers to the gene mutation that occurs in somatic cells other than sex cells. Somatic cell variation can cause changes in cell function. Specifically, the software Metcts can be used to identify the somatic cell variation of the sample.

[0049] Preferably, the optional steps include filtering the sample sequencing data to remove low-quality reads and reads containing adapter or sequencing primer sequences; and performing base quality correction on the comparison results.

[0050] S03, sample variation filtering

[0051] Specifically, variants with poor credibility are filtered out first, such as variants with fewer supporting reads and variants that appear in normal control samples at the same time; then, variants located in intergenic regions and introns and synonymous variants are filtered out based on the annotation results. Such variants will not cause changes in the protein sequence.

[0052] S04, translation of variant sequences into proteins

[0053] Specifically, among the variation results that can produce protein sequence changes after filtering, mutant transcripts are constructed based on the genome mutation information and annotation information and translated according to the translation rules to obtain mutant protein sequences.

[0054] S05, protein fragment segmentation and filtering of peptides in the normal human proteome

[0055] Specifically, the protein sequence obtained in the previous step is first divided into k-mers of smaller length. A k-mer is a set of all possible substrings of length k contained in a string. For an input protein sequence, starting from the first amino acid residue, a sliding window with a step size of 1 is used to extract sequences of fixed length k in sequence. These sequences are k-mers. More specifically, the protein sequence obtained in S04 is divided into 9 to 12 k-mers using software.

[0056] Secondly, we filter the peptides in the normal human proteome. Specifically, UniProt is a comprehensive non-redundant database that contains all public human protein sequences, and removes Kmer fragments in UniProt proteins.

[0057] S06, predicting protein fragment structures

[0058] Generally speaking, antigens first bind to MHC proteins on APCs, and then APCs present MHC-antigens to TCRs on the surface of T cells to generate immune responses. Currently, the method of screening immunologically active peptides through sequencing and bioinformatics is usually predicted by the affinity between peptides and MHC molecules. There is increasing evidence that the structural features of peptides affect T cell recognition independently of the affinity between peptides and MHC molecules. The spatial folding form of proteins is the basic unit of their biological activity and function, and can provide information about protein functions and binding characteristics. The algorithm can accurately predict the structure of peptides. Specifically, the segmented protein sequence is input using the software pep-fold, and according to the corresponding algorithm, the tertiary structure diagram of the polypeptide and the pdb file (protein three-dimensional structure data file) can be obtained.

[0059] S07, calculation of corresponding peptide structural indexes

[0060] Specifically, the structural indicators related to the polypeptide are calculated, the obtained polypeptide pdb file is imported into the Gromacs software, and the commands of the desired related structural indicators are input, and finally the result file of the corresponding structural indicators can be obtained. Through analysis, the immunogenic polypeptides can be screened out.

[0061] S08, Screening for new antigens

[0062] Usually, the mutant proteins expressed by tumor cells cannot be expressed by normal cells. These mutant proteins are processed into short peptides by proteases in cells, and then bind to human leukocyte antigens, present on the cell surface and recognized by T cells as foreign antigens. Immunogenic peptides can be screened out through the structural indicators of peptides.

[0063] First, RMSF was used to preliminarily screen the peptides. According to the research results, the RMSF value of the C-terminus of the mutant polypeptide chain with immune activity is usually less than The RMSF values ​​obtained by the above method are sorted from small to large RMSF. Generally speaking, the larger the C-terminal RMSF value, the more unstable the C-terminal of the peptide is, and the lower the efficiency of binding to the MHC protein. For this reason, RMSF can be used to preliminarily screen out peptides that may be immunogenic.

[0064] Secondly, SASA is used for further screening. Studies have shown that immunogenic peptides usually have larger SASA values ​​than non-immunogenic peptides. At the same time, immunogenic mutant peptides often have increased SASA values ​​compared to their corresponding wild-type peptides. Usually, the difference in SASA is greater than 0.2nm 2 , less than 1nm2 When the SASA difference is positively correlated with the immunogenicity of the peptide, that is, the larger the SASA difference, the stronger the immunogenicity of the mutant peptide. To this end, the SASA of the wild-type and mutant peptides were calculated respectively, and the SASA value of the mutant peptide was subtracted from the corresponding wild-type peptide SASA value, and then the obtained difference was used to sort, and the difference between 0.2 and 1 nm was selected. 2 The larger the difference between the two peptides, the greater the possibility that the corresponding mutant peptide has stronger immunogenicity.

[0065] The specific parameters provided by the software are as follows:

[0066] Use BWAMEM to align the sequencing data. The example command is:

[0067]

[0068] Among them, -R represents the alignment result header file, -t represents the number of running threads, -M represents the index file used, reference.fa represents the reference genome fasta file, and in.1.fq and in.2.fq represent sequencing data.

[0069] Use Picard MarkDuplicates to remove PCR duplicates in the alignment data. The example command is:

[0070] |6.java-jar picard.jar MarkDuplicates I=in.bam O=out.bam M=picard1.txt

[0071] Where I represents the input alignment file, O represents the output alignment file, and M represents the output result statistics file. Use GATK BaseRecalibrator to perform base quality correction, and the example command is:

[0072]

[0073] Among them, -R represents the reference genome file, -I represents the input BAM file, -O represents the output statistical result file, and --known-sites represents the known mutation file.

[0074] Use Mutect2 to identify variants. The example command is:

[0075]

[0076] Where -R represents the reference genome fasta file, -I represents the input alignment file, -tumor / -normal represents the name of the tumor / normal sample in the alignment file, and -O represents the output mutation file.

[0077] Use GATK FilterMutectCalls to perform preliminary filtering on the mutation results. The example command is:

[0078]

[0079] Where -V represents the input mutation file, and -O represents the output mutation file with the FILTER tag. Use VEP to annotate mutations. The example command is:

[0080] |15.perl vep.pl-i in.vcf-o out.txt--assembly assembly--fork 10

[0081] Where -i represents the input mutation file, -o represents the output result file, --assembly represents the reference genome version, and --fork represents the number of threads

[0082] Use Gromacs software to obtain the topology file, which contains all force field parameters (based on the selected force field) and calculate the peptide structural indexes. The example command is:

[0083] |16.gmx pdb2gmx-f model1.pdb-water tip3p

[0084] Where -f specifies the protein structure file to be processed, -water specifies the water model, and the present invention uses the TIP3P water model. This step will generate the structure file gro, the topology file top, and the position constraint file itp. Use the editconf command to create a periodic simulation box, and its example command is:

[0085] |17.gmx editconf-f conf.gro–bt dodecahedron–d 1.0–o box.gro

[0086] Among them, -f represents the input protein structure, -bt creates a rhombic dodecahedron box, because this box is close to a sphere and has the highest computational efficiency, -d sets the minimum distance from the molecule to the edge of the box, in nm, which determines the size of the box. Theoretically, in most systems, -d cannot be less than 0.9nm, -o: outputs a structure file with simulation box information.

[0087] Use gmx solvate to fill a periodic box of a given size / type with the appropriate number of solvent molecules and perform energy minimization. An example command is:

[0088]

[0089] Among them, -cp specifies the system that needs to be filled with water molecules, with a simulated protein box, -cs specifies the use of the SPC water model for filling, spc216 is the unified three-point water molecule structure of GROMACS, -p modifies the topology file of the system and adds the physical parameters of the corresponding water molecules, and -o specifies the output file after filling with water molecules.

[0090] gmx grompp (gromacs preprocessor) reads the molecular topology file, checks the validity of the file, and expands the topology from molecular description to atomic description. -f indicates the input MD parameter file, -p indicates the input topology file, -c indicates the input structure file, -r indicates the input structure file, and -o specifies the tpr file.

[0091] The gmx genion command will replace some water molecules with ions. -s indicates the input tpr file, -o indicates the input structure file, -p indicates the input topology file, -pname indicates the cation name, -np indicates the number of cations, -nname indicates the anion name, and -nn indicates the number of anions.

[0092] If molecular dynamics is started immediately, the added hydrogens in water and the broken hydrogen bond network will lead to considerable forces and structural deformations. In order to eliminate these forces, it is necessary to first run a short energy minimization. Gromacs uses a separate preprocessor, grompp, to collect the parameters, topology and coordinates into a single run input file (em.tpr), from which the simulation is then started. gmx mdrun is the main computational chemistry engine in GROMACS and can be used to calculate energy minimizations. -deffnm sets the default file name for all file options.

[0093] To avoid unnecessary protein distortions, an equilibration run is performed before the molecular dynamics simulation begins. All heavy protein atoms are constrained to their starting positions (using the file posre.itp generated earlier) and the water molecules around the protein are equilibrated. An example command is:

[0094]

[0095] After completing the two stages of equilibration, the system is now at the desired temperature and pressure. We can now release the position constraints and run the finished MD to collect data. An example command would be:

[0096] |26.gmx grompp-f run.mdp-p topol.top-c pr.gro-r pr.gro-o run.tpr

[0097] |27.gmx mdrun-v-deffnm run

[0098] The mdp file used in the above steps is an annotation file, which contains annotations for each item in the simulation process. For example, specify the integration algorithm, integration step, charge calculation method, specify the system temperature, pressure and other information. Each mdp file is adjusted accordingly according to the role of each step.

[0099] Related data show that the stability of the C-terminus of a peptide is related to whether the peptide can induce an immune response. RMSF refers to the root mean square fluctuation of the atomic position, which can be used to evaluate the stability of the C-terminus of a peptide. The example command for calculating the RMSF value of the C-terminus of a peptide-related structural indicator using gmx rmsf is:

[0100]

[0101] Where -s indicates input tpr file, -f indicates track file, -o indicates output xvgr / xmgr file, and -oq indicates protein database file.

[0102] SASA refers to the solvent accessible surface area, which is an important factor affecting protein folding. Generally speaking, the SASA of immunogenic peptides is greater than that of non-immunogenic peptides. SASA of peptide molecules can be calculated using gmx sasa. -s indicates the input structure file tpr, -f indicates the input trajectory or single configuration file, -o indicates that SASA is defined as a time function input xvg file, -odg indicates the estimation of the change of solvent free energy over time, -or calculates the average area of ​​each residue, and -oa calculates the average area of ​​each atom.

[0103] The preferred specific implementation modes and embodiments of the present invention are described in detail above, but the present invention is not limited to the above implementation modes and embodiments, and various changes can be made within the knowledge scope of those skilled in the art without departing from the concept of the present invention.

Claims

1. A method for predicting new antigens based on polypeptide structural indicators, characterized in that: The following steps are involved: S01, obtaining sample sequencing data; S02, identification of sample somatic variants; S03, sample variation filtering; S04, translation of variant sequence into protein; S05, protein fragment segmentation and filtering of peptides in the normal human proteome; S06, predicting protein fragment structure; S07, calculating the structural index of the corresponding polypeptide; S08, screening for new antigens; S08 includes the following steps: S81, preliminary screening of peptides based on RMSF; S82, further screening based on SASA; In S01, DNA sequencing is performed on the tumor sample and the normal control sample to obtain sample sequencing data, wherein the sample sequencing data includes a plurality of overlapping or partially overlapping short read sequences; in S01, no less than 500X deep sequencing data of the tumor sample and 200X deep sequencing data of the normal control are obtained; In S02, the following steps are included: S21, aligning the sequencing data of the tumor sample and the normal control sample to the reference genome, and locating the position of the short read sequence on the reference genome; S22, removing PCR duplication from the result file; S23, identifying somatic cell variation of the sample; in S02, the somatic cell variation refers to gene mutation occurring in somatic cells other than sex cells; S03 includes the following steps: S31, filtering out variants with poor credibility; S32, filtering variants located in intergenic regions and introns as well as synonymous variants based on the annotation results; In S04, among the variant results that can produce protein sequence changes after filtering, a mutant transcript is constructed according to the genome mutation information and annotation information and translated according to the translation rules to obtain a mutant protein sequence; In S05, the obtained protein sequence is segmented into Kmers with a length of 9 to 12, and the peptides in the normal human proteome are filtered; In S06, the segmented protein sequence is input into the software pep-fold, and according to the corresponding algorithm, the tertiary structure diagram of the polypeptide and the protein three-dimensional structure data file are obtained; In S07, the result files of structural indicators were obtained using Gromacs software, and peptides with immunogenicity were screened out through analysis; The difference of screening SASA is between 0.2 and 1 nm 2 between the peptides.

2. The method according to claim 1, characterized in that In S01, whole exome sequencing or panel capture sequencing is used.

3. The method according to claim 1, characterized in that It also includes filtering the sample sequencing data, removing low-quality reads and reads containing adapter or sequencing primer sequences, and performing base quality correction on the comparison results.

4. The method according to claim 1, characterized in that: The variants with poor credibility include variants with fewer supporting reads and variants that also appear in normal control samples.

Images