Method for rapidly screening human cell exosome key marker based on RNA-seq

Abstract

The invention discloses a method for rapidly screening key markers of human cell exosomes based on RNA-seq, which comprises the following steps: S1, respectively carrying out RNA-seq sequencing on a normal cell line, a tumor cell line varied from the normal cell line and respective cell exosomes of the normal cell line and the tumor cell line to obtain a target expression matrix; s2, respectively screening common genes of a normal cell line and a cell exosome of the normal cell line and common genes of a tumor cell line and a cell exosome of the tumor cell line by using a computer script to obtain a new expression profile matrix; s3, carrying out weighted gene co-expression network analysis on the new expression spectrum matrix, and identifying key genes; s4, extracting membrane protein verified by experiments from the identified key genes; and S5, performing significant differential expression analysis, and extracting significant differential expression genes. According to the invention, exosome analysis can be rapidly carried out on a certain type of pathological cells, and cell exosome key marking is completed.

Description

Technical Field

[0001] The present invention relates to basic medical research and provides a method for rapidly screening key markers of human cell exosomes based on RNA-seq. Background Art

[0002] Exosomes, one of the three major markers of liquid biopsy, have been gaining increasing attention in recent years. Whether in the field of oncology, neuroscience or metabolism, a large number of exosome-related articles have been published to explore the functions and mechanisms of exosomes.

[0003] Exosomes are a type of extracellular vesicles with a size of 30-150nm, saucer-shaped, and wrapped by a lipid bilayer membrane. They are widely present in cell culture supernatants and various body fluids, including blood, lymph, saliva, urine, follicular fluid, semen, pericardial fluid, milk, etc. Exosomes contain a variety of proteins (such as TSG101, HSP70, etc.), nucleic acids (DNA, RNA) and lipids (cholesterol, sphingomyelin, etc.). These contents vary greatly depending on the cells that release the exosomes. It can be seen that exosomes are extracellular vesicles produced by all cells, carrying nucleic acids, proteins, lipids and metabolites. They are the medium for close communication between healthy and diseased cells, affecting all aspects of cell biology and are of great research value.

[0004] In order to study the mechanism and function of exosomes, it is first necessary to separate exosomes from the sample. The most recognized separation method is ultracentrifugation, followed by identification by transmission electron microscopy (TEM), nanoparticle tracking analysis (NTA), Western Blot, etc. After obtaining the exosomes, the exosomes are labeled, the uptake of the exosomes is observed, and the functions of the exosomes are further explored, such as tissue repair, tumor metastasis, metabolic reconstruction, etc. According to the current mainstream view, exosome labeling is an indispensable step in exploring the function of exosomes. First, it is necessary to observe that the exosomes are absorbed by the receptor cells or reach the target tissues and organs, and then detect the changes in the corresponding indicators in the receptor cells, tissues or organs to see whether the exosomes have certain functions. After confirming that the exosomes have functions, high-throughput sequencing or gene chips are used to mine differentially expressed molecules and further study the mechanism.

[0005] The current methods for labeling exosomes are mainly: 1. Lipophilic dye labeling: In the currently published exosome articles, most exosomes are labeled with lipophilic dyes, which are widely used both in vivo and in vitro. Lipophilic dyes are mainly divided into two categories. The first category is PKH67 (green fluorescence) / PKH26 (red fluorescence). Since they can stably bind to the lipid bilayer membrane of exosomes, they have better staining effects and are widely used. The second category is the Di series of lipophilic dyes, including DiI (orange fluorescence), DiO (green fluorescence), DiD (red fluorescence), and DiR (deep red fluorescence). Among them, the infrared fluorescence of DiR can penetrate cells and tissues and is used for tracing in in vivo imaging. 2. Lentivirus-mediated CD63-GFP expression: The expression elements of the specific protein CD63 and green fluorescent protein GFP of the exosome are constructed into plasmids and then packaged into lentiviruses. The cells are then infected with this lentivirus, so that the exosomes secreted by the cells have green fluorescence. 3. Labeling with gold nanoparticles (GNPs): Intranasal administration of exosomes labeled with gold nanoparticles can not only track exosomes non-invasively, but also promote the accumulation of exosomes in the brain and prolong their existence in the lesion area, enhancing the effect of in vivo neuroimaging. Therefore, this exosome labeling technology can be used as a powerful diagnostic tool for a variety of brain diseases and may enhance exosome-based neurorecovery therapy. 4. Labeling with near-infrared (NIR) fluorescent probes: This is a method for endogenous detection of NIR fluorescent-labeled exosomes using amphiphilic probes, without the need for immunolabeling or synthesis or chromatographic treatment of exosomes. Since fluorescent probes with NIR wavelengths can penetrate tissues in the body, they provide the benefit of direct conversion from cells to in vivo use. 5. Labeling with Gaussialuciferase (gluc) and truncated lactadherin: gLuc-lactadherin is composed of a fusion expression plasmid, and this plasmid is transfected into cells, so that the exosomes secreted by the cells carry strong luciferase activity. The labeled exosomes were injected intravenously into mice to observe the distribution of exosomes in tissues.

[0006] RNA-Seq sequencing (transcriptome sequencing) is a transcriptomics research method based on the second-generation sequencing technology: first, all the transcribed RNA of the biological sample is extracted and mRNA is enriched, then reverse transcribed into cDNA for high-throughput sequencing, and then the fragments are spliced ​​and assembled on this basis to obtain individual transcripts, thereby forming a global understanding of the gene expression status of the current developmental state of the biological sample. Summary of the invention

[0007] Purpose of the invention: The purpose of the present invention is to provide an RNA-seq-based exosome research method that is different from the existing research route to quickly screen key markers of human cell exosomes.

[0008] Technical solution: The method for rapid screening of key markers of human cell exosomes based on RNA-seq of the present invention comprises the following steps:

[0009] S1. Perform RNA-seq sequencing on normal cell lines, tumor cell lines mutated from normal cell lines, and their respective cell exosomes, and perform data preprocessing on the sequencing results to obtain a target expression matrix;

[0010] S2. Using computer scripts, we screened the common genes of exosomes between normal cell lines and normal cell lines, and the common genes of exosomes between tumor cell lines and tumor cell lines, respectively, to obtain expression profile data of two groups of exosomes, and merged the two groups of data to obtain a new expression profile matrix;

[0011] S3, constructing a WGCNA network, reading the new expression spectrum matrix obtained in step S2 to perform weighted gene co-expression network analysis and identify key genes;

[0012] S4. Extract experimentally verified membrane proteins from the identified key genes;

[0013] S5. Use DESeq2 to perform significant differential expression analysis, extract significantly differentially expressed genes, and draw a volcano plot of significantly differentially expressed genes in the exosomes of normal cell lines and tumor cell lines.

[0014] A further preferred technical solution of the present invention is that when there are multiple types of cells associated with the same tumor disease in step S1, the different cells are divided into multiple groups of cell lines, the corresponding tumor cell lines, and the respective cell exosomes are also divided into multiple groups, and RNA-seq sequencing is performed separately.

[0015] Further preferably, in step S2, an R script is used to call a function to screen common genes.

[0016] Further preferably, in step S4, firstly, a database of exosome-related markers is downloaded, and then the database is compared with the key genes identified in step S3 to extract membrane proteins that have been experimentally verified.

[0017] Preferably, after obtaining the significantly differentially expressed genes, gene function enrichment analysis is performed by enrichment.

[0018] Beneficial effects: The present invention provides a method for rapidly screening key markers of human cell exosomes based on RNA-seq, which can rapidly perform exosome analysis on a certain type of diseased cells without the need for biomarkers, and directly obtain significantly differentially expressed genes in the exosomes of tumor cells to complete key markers of exosomes; on this basis, significantly differentially expressed genes can be enriched by enrichment to perform functional analysis. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] Figure 1 The volcano map of significantly differentially expressed genes in exosomes drawn in step S5 of the embodiment of the present invention;

[0020] Figure 2 This is a functional enrichment analysis diagram of step S6 in an embodiment of the present invention. DETAILED DESCRIPTION

[0021] The technical solution of the present invention is described in detail below with reference to the accompanying drawings, but the protection scope of the present invention is not limited to the embodiments.

[0022] Embodiment: This embodiment provides a specific application of a method for rapid screening of key markers of human cell exosomes based on RNA-seq. In this embodiment, glioblastoma is taken as an example. The cells related to glioblastoma are: human brain glial cells and astrocytes, so this embodiment is divided into 4 groups of cell lines. Each group is specifically as shown in the following table. HEB is a human brain glial cell, HA1800 is a normal astrocyte, HEBexo is a human brain glial cell exosome, HA1800exo is a normal astrocyte cell exosome, T98G and U251 are human glioma cells, T98Gexo and U251exo are human glioma cell exosomes.

[0023] Table 1: Glioblastoma sample grouping table

[0024]

[0025]

[0026] The specific steps are as follows:

[0027] S1. Perform RNA-seq sequencing on each group of samples in Table 1, and perform data preprocessing on the sequencing results to obtain the target expression matrix;

[0028] After upstream RNA-seq preprocessing, readcounts, Allcount.txt, are obtained, which are separated by tabs. After annotation, all genes are obtained. The first row symbol, HEB_1, HEB_2, HEB_3, HA1800_1, HA1800_2, HA1800_3, HEBexo_1, HEBexo_2, HEBexo_3, HA1800exo_1, HA1800exo_2, HA1800exo_3, T98G_1, T98G_2, T98G_3, U251_1, U251_2, U251_3, T98G_1exo, T98G_2exo, T98G_3exo, U251exo_1, U251exo_2, U251exo_3 are sample names. Each line in the first column is a gene name, which exists in symbol.

[0029] #Delete all rows except symbol where all other rows are 0

[0030] awk'{for(i=2;i<=NF;i++)if($i!=0){print;next}}'Allcount.txt>Allcount_filtered.txt

[0031] explain:

[0032] The awk command is used to process text files, where NF represents the number of fields in the current line, and $i! = 0 means that the i-th field is not equal to 0.

[0033] # Sort the data except the first row by the gene name in the first column

[0034] tail-n+2Allcount_filtered.txt|sort-k1,1>Allcount_filtered_sorted.txt

[0035] #Merge sort results and table headers

[0036] (cat<(head-n 1Allcount_filtered.txt)Allcount_filtered_sorted.txt)>Allcount_filtered_sorted_with_header.txt

[0037] # Delete duplicate genes

[0038] awk-F'\t"NR==1{print$0}NR>1{if(!($1in genes)){genes[$1]=$0;print$0}}'Allcount_filtered_sorted_with_header.txt>Allcount_filtered_unique.txt

[0039] explain:

[0040] tail-n+2Allcount_filtered.txt means output all rows except the first row

[0041] sort-k1,1 means sorting by gene name in the first column

[0042] >Allcount_filtered_sorted.txt means saving the sorting results to a file named Allcount_filtered_sorted.txt

[0043] (cat<(head-n 1Allcount_filtered.txt)Allcount_filtered_sorted.txt) means merging the header and sorting results

[0044] >Allcount_filtered_sorted_with_header.txt means saving the merged results to a file named Allcount_filtered_sorted_with_header.txt

[0045] -F'\t' means use tab as field separator

[0046] NR == 1{print$0} means outputting the first line directly, i.e. the table header

[0047] NR>1{if(!($1in genes)){genes[$1]=$0;print$0}} means that for each row of data, if the gene name in the first column of the row is not in the genes array, the row of data is added to the genes array and the row of data is output. Otherwise, no operation is performed. In this way, the first occurrence of each gene name in the table can be retained, and the subsequent rows of duplicate gene names can be deleted.

[0048] The output is finally saved to a file called Allcount_filtered_unique.txt

[0049] In the final generated Allcount_filtered_unique.txt file, the first row remains unchanged, and the subsequent rows are sorted according to the gene names in the first column, and the rows with duplicate gene names are deleted.

[0050] S2. Using computer scripts, we screened the common genes of exosomes between normal cell lines and normal cell lines, and the common genes of exosomes between tumor cell lines and tumor cell lines, respectively, to obtain expression profile data of two groups of exosomes, and merged the two groups of data to obtain a new expression profile matrix;

[0051] #Use R script to screen the transcriptome sequencing data of glial cells HEB and HA1800 (6 samples) and cell exosomes HEB and HA1800, obtain common genes, and extract the expression profile data of cell exosomes HEB and HA1800A.

[0052] The transcriptome data of T98G_1, T98G_2, T98G_3, U251_2, U251_3 and cell exosomes U251exo_1, U251exo_2, U251exo_3 were screened to obtain common genes, and the expression profile data B of cell exosomes of glioma cell lines were extracted. A new expression profile matrix was obtained, with the row name as gene symbol and the column name as cell exosome samples HEBexo_1, HEBexo_2, HEBexo_3, HA1800exo_1, HA1800exo_2, HA1800exo_3, T98G_1exo, T98G_2exo, T98G_3exo, U251exo_1, U251exo_2, U251exo_3.

[0053] #Read expression matrix

[0054] expr_matrix<-read.table("Allcount_filtered_unique.txt", header=TRUE, row.names=1)

[0055] #View the dimensions of the expression matrix

[0056] dim(expr_matrix)

[0057] #View the first few rows of the expression matrix

[0058] head(expr_matrix)

[0059] Among them, the read.table() function is a function for reading text files, header = TRUE means that the first row is the column name, and row.names = 1 means that the first column is the row name. After reading, you can use the dim() function to view the dimension of the expression matrix, and use the head() function to view the first few rows of the expression matrix.

[0060] #Read expression matrix

[0061] expr_matrix<-read.table("Allcount_filtered_unique.txt", header=TRUE, row.names=1)

[0062] #Filter out the sample list of glial cells and exosomes

[0063] HEB_cols<-grep("^HEB|HEBexo",colnames(expr_matrix))

[0064] HA1800_cols<-grep("^HA1800|HA1800exo",colnames(expr_matrix))

[0065] #Extract common gene symbols

[0066] common_genes<-rownames(expr_matrix)

[0067] #Extract the expression profile data of HEB and HA1800 from cell exosomes A

[0068] A<-expr_matrix[common_genes,c(HEB_cols,HA1800_cols)]

[0069] # Check the dimension of A

[0070] dim(A)

[0071] # View the first few rows of A

[0072] head(A)

[0073] In this script, we first use the read.table() function to read the expression matrix data, and then use the grep() function to filter out the sample column names of glial cells and cell exosomes. Then we extract the gene names of all rows as common gene symbols. Then we use the common_genes variable to extract the expression profile data A of cell exosomes HEB and HA1800. Finally, we use the dim() function to view the dimensions of A, and the head() function to view the first few rows of A.

[0074] #Filter out the sample column names of T98G and U251

[0075] T98G_cols<-grep("^T98G|T98G_\\d+exo",colnames(expr_matrix))

[0076] U251_cols<-grep("^U251|U251_\\d+exo",colnames(expr_matrix))

[0077] #Extract common gene symbols

[0078] common_genes<-intersect(common_genes,rownames(expr_matrix[,c(T98G_cols,U251_cols)]))

[0079] #Extract the expression profile data of exosomes from glioma cell lines B

[0080] B<-expr_matrix[common_genes,c(T98G_cols,U251_cols)]

[0081] # Check the dimension of B

[0082] dim(B)

[0083] # View the first few lines of B

[0084] head(B)

[0085] First, use the grep() function to filter out the sample column names of T98G and U251. Then use the intersect() function to extract the intersection with the previously filtered common gene symbol as a new common gene symbol. Finally, use the new common gene symbol to extract the expression profile data B of the cell exosomes of the glioma cell line. Finally, use the dim() function to view the dimension of B, and use the head() function to view the first few rows of B.

[0086] #Extract the gene symbols of A and B

[0087] genes_A<-rownames(A)

[0088] genes_B<-rownames(B)

[0089] # Take the intersection of the gene symbols of A and B

[0090] common_genes<-intersect(genes_A,genes_B)

[0091] #Re-extract the expression data of A and B, and only keep the intersection

[0092] A<-A[common_genes,]

[0093] B<-B[common_genes,]

[0094] #Merge A and B into a new expression matrix

[0095] new_expr_matrix<-cbind(A,B)

[0096] #Set new column names

[0097] colnames(new_expr_matrix)<-c(paste0("HEBexo_",1:3),paste0("HA1800exo_",1:3),paste0("T98G_",1:3,"exo"),paste0("U251_",1:3,"exo"))

[0098] #View the dimensions of the new expression matrix

[0099] dim(new_expr_matrix)

[0100] #View the first few rows of the new expression matrix

[0101] head(new_expr_matrix)

[0102] In this script, first use the rownames() function to extract the gene symbols of A and B, and use the intersect() function to get the intersection. Next, re-extract the expression data of A and B, and keep only the intersection. Then use the cbind() function to merge A and B into a new expression profile matrix. Use the paste0() function to set the new column names to cell exosome samples HEBexo_1, HEBexo_2, HEBexo_3, HA1800exo_1, HA1800exo_2, HA1800exo_3, T98G_1exo, T98G_2exo, T98G_3exo, U251exo_1, U251exo_2, U251exo_3. Finally, use the dim() function to view the dimensions of the new expression matrix, and use the head() function to view the first few rows of the new expression profile matrix.

[0103] #Save the new expression matrix as a file

[0104] write.table(new_expr_matrix,file="new_expr_matrix.txt",sep="\t",quote=FALSE,col.names=NA)

[0105] In this script, the new expression matrix is ​​saved as a file using the write.table() function. The first argument is the data frame to be saved, the second argument is the file name, the third argument is the delimiter (here we use a tab), the fourth argument is whether quotes should be added around the values, and the fifth argument is whether the column names should be output in the first row.

[0106] ###Extract gene length and convert count to tpm

[0107] file<-"gencode.v22.annotation.gtf"

[0108] eff_length<-get_eff_len(file=file,

[0109] #Use mode 2 to calculate the total length of all exons of all transcripts

[0110] mode = 'gene',

[0111] #Use mode 1 to calculate the longest transcript length

[0112] #mode=='transcript',

[0113] feature="exon",

[0114] gene_id = 'gene_name',

[0115] transcript_id="transcript_id")

[0116] eff_length1<-data.frame(eff_length)

[0117] write.table(eff_length1,"eff_length1.csv",sep=",",quote=F,row.names=TRUE)

[0118] S3, constructing a WGCNA network, reading the new expression spectrum matrix obtained in step S2 to perform weighted gene co-expression network analysis and identify key genes;

[0119] #Read new expression matrix file

[0120] new_expr_matrix<-read.table("new_expr_matrix.txt", header=TRUE, row.names=1, sep="\t")

[0121] #Install and load the WGCNA package

[0122] install.packages("WGCNA")

[0123] library(WGCNA)

[0124] #Convert data format

[0125] new_expr_matrix<-t(new_expr_matrix)

[0126] #Build WGCNA network

[0127] powers<-c(1:10)

[0128] sft<-pickSoftThreshold(new_expr_matrix,powerVector=powers,verbose=5)

[0129] power <- sft$powerEstimate

[0130] net<-blockwiseModules(new_expr_matrix,power=power,TOMType="unsigned",mergeCutHeight=0.25,minModuleSize=30,

[0131] reassignThreshold=0,networkType="signed")

[0132] #Get the key genes of the modules where T98Gexo and U251exo are located

[0133] T98Gexo_module<-names(which(net$colors=="turquoise")%in%grep("T98Gexo_",colnames(new_expr_matrix)))

[0134] U251exo_module<-names(which(net$colors=="turquoise")%in%grep("U251exo_",colnames(new_expr_matrix)))

[0135] key_genes<-c(T98Gexo_module,U251exo_module)

[0136] In this script, first use the read.table() function to read the new expression matrix file. Next, use the install.packages() function to install the WGCNA package, and use the library() function to load the WGCNA package. Then, use the t() function to convert the new expression matrix into a format with samples as rows and genes as columns. Next, use the pickSoftThreshold() function to select the best power value, and use the blockwiseModules() function to build the WGCNA network. Then, according to the color of the module where the glioma cell exosomes T98Gexo and U251exo are located (here assumed to be turquoise), use the grep() function to extract the genes of these two modules from the column names, and merge the results into a vector. Finally, print the key genes.

[0137] # Output key genes

[0138] key_genes

[0139] #Save key genes to file

[0140] write.table(key_genes,file="key_genes.txt",row.names=FALSE,col.names=FALSE,quote=FALSE)

[0141] We use the write.table() function to save the key genes to a file called key_genes.txt and set options using the row.names, col.names, and quote parameters to control the format of the file.

[0142] Extract the expression matrix in new_expr_matrix.txt according to the symbol in key_genes.txt to form key_expr.txt

[0143] #Read new_expr_matrix.txt and key_genes.txt files

[0144] new_expr_matrix<-read.table("new_expr_matrix.txt", header=TRUE, row.names=1, sep="\t")

[0145] key_genes<-read.table("key_genes.txt",header=FALSE,col.names=c("symbol"))

[0146] #Match the symbol in key_genes with the row name in new_expr_matrix and extract the corresponding row

[0147] key_expr_matrix<-new_expr_matrix[match(key_genes$symbol,rownames(new_expr_matrix)),]

[0148] #Write key_expr.txt file

[0149] write.table(key_expr_matrix,file="key_expr.txt",quote=FALSE,sep="\t",col.names=NA)

[0150] In this script, we first use the read.table() function to read the new_expr_matrix.txt and key_genes.txt files. Then, we use the match() function to match the symbols in key_genes with the row names in new_expr_matrix and extract the corresponding rows to get the key_expr_matrix. Next, we use the write.table() function to save the key_expr_matrix to a file called key_expr.txt and use the sep, quote, and col.names parameters to set options to control the format of the file.

[0151] S4, download the exosome-related marker database, compare the database with the key genes identified in step S3, and extract the membrane proteins that have been experimentally verified;

[0152] Download EXOCARTA_PROTEIN_MRNA_DETAILS_5.txt from the official website

[0153] #Read the EXOCARTA_PROTEIN_MRNA_DETAILS_5.txt file

[0154] data<-read.table("EXOCARTA_PROTEIN_MRNA_DETAILS_5.txt",sep="\t",header=T,check.names=F,stringsAsFactors=F,row.names=1)

[0155] #Extract rows where SPECIES is Homo_sapiens

[0156] filtered_data<-subset(data,SPECIES=="Homo_sapiens")

[0157] #Read gene_symbol.txt file

[0158] gene_symbols<-read.table("gene_symbol.txt",header=TRUE,sep="\t")

[0159] #Extract the lines that match the gene symbol in the gene_symbol.txt file

[0160] matched_data<-subset(filtered_data,GENE_SYMBOL%in%gene_symbols$gene_symbol)

[0161] #Determine whether it is a duplicate row based on the GENE_SYMBOL column

[0162] duplicate_rows<-duplicated(matched_data$GENE_SYMBOL)

[0163] #Keep the first occurrence of the line

[0164] matched_data_unique<-matched_data[! duplicate_rows,]

[0165] # Output the results to a new file

[0166] write.table(matched_data_unique,"matched_data_unique.txt",sep="\t",quote=FALSE,row.names=FALSE)

[0167] #Only keep the GENE_SYMBOL and EXPERIMENT_ID METHODS columns

[0168] final_data<-matched_data_unique[,c("GENE_SYMBOL","METHODS")]

[0169] # Output the results to a new file

[0170] write.table(final_data,"output.txt",sep="\t",quote=FALSE,row.names=FALSE)

[0171] library(readr)

[0172] #Read output.txt file

[0173] output_data<-read.table("output.txt",sep="\t",header=T,stringsAsFactors=F)

[0174] #Extract the GENE_SYMBOL column

[0175] gene_symbols<-output_data$GENE_SYMBOL

[0176] #Read All_reads_counts.csv file

[0177] all_reads_counts<-read.csv("All_reads_counts.csv",header=T)

[0178] #First, match the elements in gene_symbols with the Geneid column in the all_reads_counts data frame to find the matching row index. This can be achieved by using the %in% operator. For example:

[0179] rows<-all_reads_counts$Geneid%in%gene_symbols

[0180] #Use the row index extracted in the previous step to extract the matching rows from the all_reads_counts data frame. For example:

[0181] matching_rows1<-all_reads_counts[rows,]

[0182] #The variable matching_rows will contain the rows in the all_reads_counts data frame that match gene_symbols.

[0183] #Read the "exo_deseq2_test_result.exo_HEB_HA1800_vs_exo_T98G_U251.txt" file and specify the correct delimiter and column name. For example:

[0184] result_data<-read.table("exo_deseq2_test_result.exo_HEB_HA1800_vs_exo_T98G_U251.txt",sep="\t",header=TRUE)

[0185] #Please make sure that the specified file path and delimiter are correct and that the file contains column names.

[0186] #Use the elements in gene_symbols to match the id column in the result_data data frame to find the matching row index. This can be achieved by using the %in% operator. For example:

[0187] rows2<-result_data$id%in%gene_symbols

[0188] #Use the row index extracted in the previous step to extract the matching rows from the result_data data frame. For example:

[0189] matching_rows2<-result_data[rows2,]

[0190] #merge

[0191] merged_data<-merge(matching_rows1,matching_rows2,by.x=1,by.y=1)

[0192] #Output

[0193] write.csv(merged_data,file="match_expr_exoNormal_vs_exoGlioma.csv",row.names=FALSE)

[0194] S5. Use DESeq2 to perform significant differential expression analysis, extract significantly differentially expressed genes, and draw a volcano map of significantly differentially expressed genes in the exosomes of normal cell lines and tumor cell lines, such as Figure 1 As shown;

[0195] ###Significant differential expression analysis using DESeq2

[0196] The key_expr.txt file contains row names of gene symbol and column names of exosome samples HEBexo_1, HEBexo_2, HEBexo_3, HA1800exo_1, HA1800exo_2, HA1800exo_3, T98G_1exo, T98G_2exo, T98G_3exo, U251exo_1, U251exo_2, U251exo_3. HEBexo_1, HEBexo_2, HEBexo_3, HA1800exo_1, HA1800exo_2, HA1800exo_3 are control, T98G_1exo, T98G_2exo, T98G_3exo, U251exo_1, U251exo_2, U251exo_3 are glioma, which is the grouping information.

[0197] #Read the key_expr.txt file

[0198] key_expr_matrix<-read.table("key_expr.txt",header=TRUE,row.names=1,sep="\t")

[0199] #Build the experimental design

[0200] group<-factor(c(rep("control",6),rep("glioma",6)))

[0201] colData<-data.frame(sample=colnames(key_expr_matrix),group)

[0202] #Define group information

[0203] control<-c("HEBexo_1","HEBexo_2","HEBexo_3","HA1800exo_1","HA1800exo_2","HA1800exo_3")

[0204] glioma<-c("T98G_1exo","T98G_2exo","T98G_3exo","U251exo_1","U251exo_2","U251exo_3")

[0205] # Group samples

[0206] control_samples<-colData$sample[colData$group=="control"]

[0207] glioma_samples<-colData$sample[colData$group=="glioma"]

[0208] #Use DESeq2 for differential expression analysis

[0209] library(DESeq2)

[0210] dds<-DESeqDataSetFromMatrix(countData=key_expr_matrix,colData=colData,design=~group)

[0211] dds<-DESeq(dds)

[0212] res<-results(dds,contrast=c("group","glioma","control"))

[0213] #Extract significantly differentially expressed genes

[0214] sig_genes<-rownames(res)[which(res$padj<0.05&abs(res$log2FoldChange)>2)]

[0215] ##Volcano plot of significantly differentially expressed genes in exosomes of glial cells and glioma cells

[0216] library(ggplot2)

[0217] library(ggpubr)

[0218] library(ggthemes)

[0219] library(cowplot)

[0220] deg.data1<-read.csv("exo_deseq2_test_result.exo_HEB_HA1800_vs_exo_T98G_U251.csv",header=T)

[0221] #Draw p1

[0222] {

[0223] # Display the first 6 lines of the file

[0224] head(deg.data1)

[0225] # Perform log10 transformation on p after correction of differentially expressed genes

[0226] deg.data1$logP<--log10(deg.data1$adj.P.Val)

[0227] #Add a new column group

[0228] deg.data1$Group="not-significant"

[0229] #Set genes with adj.P.Val<0.05,logFC>2 as significantly upregulated genes

[0230] #Set genes with adj.P.Val<0.05,logFC<-2 as significantly upregulated genes

[0231] deg.data1$Group[which((deg.data1$adj.P.Val<0.05)&(deg.data1$logFC>=2))]="up-regulated"

[0232] deg.data1$Group[which((deg.data1$adj.P.Val<0.05)&(deg.data1$logFC<=-2))]="down-regulated"

[0233] # View the number of up-regulated and down-regulated genes

[0234] table(deg.data1$Group)

[0235] #Add a new label column

[0236] deg.data1$Label=""

[0237] #Sort the p values ​​of genes with significant differential expression from small to large

[0238] deg.data1<-deg.data1[order(deg.data1$adj.P.Val),]

[0239] #Select the 10 most highly expressed genes with the smallest q

[0240] up.genes<-head(deg.data1$Gene[which(deg.data1$Group=="up-regulated")],10)

[0241] #Select the 10 genes with the smallest q for low expression

[0242] down.genes<-head(deg.data1$Gene[which(deg.data1$Group=="down-regulated")],10)

[0243] #Merge upgene and downgene into label

[0244] deg.top10.genes<-c(as.character(up.genes),as.character(down.genes))

[0245] head(deg.top10.genes)

[0246] deg.data1$Label[match(deg.top10.genes,deg.data1$Gene)]<-deg.top10.genes

[0247] }

[0248] #Add a new label column

[0249] #deg.data$Label=""

[0250] #Sort the p values ​​of genes with significant differential expression from small to large

[0251] #deg.data<-deg.data[order(deg.data$adj.P.Val),]

[0252] #Select the 20 most highly expressed genes with the smallest q

[0253] #up.genes<-head(deg.data$Gene[which(deg.data$Group=="up-regulated")],20)

[0254] #Select the 20 genes with the smallest q for low expression

[0255] #down.genes<-head(deg.data$Gene[which(deg.data$Group=="down-regulated")],20)

[0256] #Merge upgene and downgene into label

[0257] #deg.top20.genes<-c(as.character(up.genes),as.character(down.genes))

[0258] #head(deg.top20.genes)

[0259] #deg.data$Label[match(deg.top20.genes,deg.data$Gene)]<-deg.top20.genes

[0260] #}

[0261] library(RColorBrewer)

[0262] pdf(file="exo_deseq2_volcano.pdf", width=8, height=6)

[0263] #theme_set(theme_bw())

[0264] #Change the axis labeling to be more beautiful

[0265]

[0266] Blue indicates key genes that are significantly highly expressed in glioma cell exosomes. The top 10 are MAP1A, ARRDC4, SH3PXD2A, CCND1, RAB3B, ORAOV1, TAGLN, KIAA0556, SAMD4A, and DPP4.

[0267] S6. After obtaining significantly differentially expressed genes, gene function enrichment analysis was performed through enrichment. The results of enrichment analysis are as follows: Figure 2 .

[0268] #Perform gene function enrichment analysis

[0269] library(clusterProfiler)

[0270] #GO Enrichment

[0271] enrichGO<-enrichGO(sig_genes,OrgDb="org.Hs.eg.db",keyType="SYMBOL",ont="BP",pvalueCutoff=0.05, qvalueCutoff=0.05)

[0272] #keggenrichment

[0273] enrichKEGG<-enrichKEGG(sig_genes,organism="hsa",pvalueCutoff=0.05, qvalueCutoff=0.05)

[0274] #GSEA analysis

[0275] gse<-gseGO(sig_genes,OrgDb="org.Hs.eg.db",ont="BP",pvalueCutoff=0.05,nPerm=1000)

[0276] # Output enrichment results

[0277] enrichGO

[0278] enrichKEGG

[0279] gse

[0280] #Write the significantly differentially expressed genes into a txt file

[0281] write.table(sig_genes,"sig_genes.txt",col.names=FALSE,row.names=TRUE,quote=FALSE,sep="\t")

[0282] #Write GO enrichment analysis results into a txt file

[0283] write.table(enrichGO,"enrichGO.txt",col.names=TRUE,row.names=FALSE,quote=FALSE,sep="\t")

[0284] #Write the kegg enrichment analysis results into a txt file

[0285] write.table(enrichKEGG,"enrichKEGG.txt",col.names=TRUE,row.names=FALSE,quote=FALSE,sep="\t")

[0286] #Write GSEA analysis results into a txt file

[0287] write.table(gse,"gse.txt",col.names=TRUE,row.names=FALSE,quote=FALSE,sep="\t")

[0288] As described above, although the present invention has been shown and described with reference to specific preferred embodiments, it should not be construed as limiting the present invention itself. Various changes in form and details may be made without departing from the spirit and scope of the present invention as defined in the appended claims.

Claims

1. A method for rapid screening of key markers of human cell exosomes based on RNA-seq, characterized in that: The steps include: S1. Perform RNA-seq sequencing on normal cell lines, tumor cell lines mutated from normal cell lines, and their respective cell exosomes, and perform data preprocessing on the sequencing results to obtain a target expression matrix; S2. Using computer scripts, we screened the common genes of exosomes between normal cell lines and normal cell lines, and the common genes of exosomes between tumor cell lines and tumor cell lines, respectively, to obtain expression profile data of two groups of exosomes, and merged the two groups of data to obtain a new expression profile matrix; S3, constructing a WGCNA network, reading the new expression spectrum matrix obtained in step S2 to perform weighted gene co-expression network analysis and identify key genes; S4. Extract experimentally verified membrane proteins from the identified key genes; S5. Use DESeq2 to perform significant differential expression analysis, extract significantly differentially expressed genes, and draw a volcano plot of significantly differentially expressed genes in the exosomes of normal cell lines and tumor cell lines.

2. The method for rapid screening of key markers of human cell exosomes based on RNA-seq according to claim 1, characterized in that: When there are multiple types of cells associated with the same tumor disease in step S1, the different cells are divided into multiple groups of cell lines, and the corresponding tumor cell lines and their respective cell exosomes are also divided into multiple groups, and RNA-seq sequencing is performed separately.

3. The method for rapid screening of key markers of human cell exosomes based on RNA-seq according to claim 1, characterized in that: In step S2, the R script is used to call the function to screen the common genes.

4. The method for rapid screening of key markers of human exosomes based on RNA-seq according to claim 1, characterized in that: In step S4, first, the exosome-related marker database is downloaded, and then the database is compared with the key genes identified in step S3 to extract membrane proteins that have been experimentally verified.

5. The method for rapid screening of key markers of human cell exosomes based on RNA-seq according to claim 1, characterized in that: After obtaining the significantly differentially expressed genes, gene function enrichment analysis was performed through enrichment.

Images