Similarity-calculation-by-sequence-alignment.R

#Similarity calculation by sequence alignment
rm(list = ls())
###############################input data 
dir_path <- "C:\\Users\\liyix\\OneDrive\\Desktop\\2024\\"
dir_path_name <- list.files(pattern = ".*",dir_path,full.names = T, recursive = T)
dir_path_name

## [1] "C:\\Users\\liyix\\OneDrive\\Desktop\\2024\\Human_kinase_domain.fasta" 
## [2] "C:\\Users\\liyix\\OneDrive\\Desktop\\2024\\Human_kinase_protein.fasta"

library("protr")

## Warning: 程辑包'protr'是用R版本4.2.3 来建造的

#Data Sample
domain <- readFASTA(grep("Human_kinase_domain.fasta",dir_path_name,value = T))[1]
domain

## $TTBK2_Hsap
## [1] "WKVLRKIGGGGFGEIYDALDMLTRENVALKVESAQQPKQVLKMEVAVLKKLQGKDHVCRFIGCGRNDRFNYVVMQLQGRNLADLRRSQSRGTFTISTTLRLGRQILESIESIHSVGFLHRDIKPSNFAMGRFPSTCRKCYMLDFGLARQFTNSCGDVRPPRAVAGFRGTVRYASINAHRNREMGRHDDLWSLFYMLVEFVVGQLPWRKIKDKEQVGSIKERYDHRLMLKHLPPEFSIFLDHISSLDYFTKPDYQLL"

nchar(domain)

## TTBK2_Hsap 
##        256

protein <- readFASTA(grep("Human_kinase_protein.fasta",dir_path_name,value = T))[1]
protein

## $TTBK2_Hsap
## [1] "MSGGGEQLDILSVGILVKERWKVLRKIGGGGFGEIYDALDMLTRENVALKVESAQQPKQVLKMEVAVLKKLQGKDHVCRFIGCGRNDRFNYVVMQLQGRNLADLRRSQSRGTFTISTTLRLGRQILESIESIHSVGFLHRDIKPSNFAMGRFPSTCRKCYMLDFGLARQFTNSCGDVRPPRAVAGFRGTVRYASINAHRNREMGRHDDLWSLFYMLVEFVVGQLPWRKIKDKEQVGSIKERYDHRLMLKHLPPEFSIFLDHISSLDYFTKPDYQLLTSVFDNSIKTFGVIESDPFDWEKTGNDGSLTTTTTSTTPQLHTRLTPAAIGIANATPIPGDLLRENTDEVFPDEQLSDGENGIPVGVSPDKLPGSLGHPRPQEKDVWEEMDANKNKIKLGICKAATEEENSHGQANGLLNAPSLGSPIRVRSEITQPDRDIPLVRKLRSIHSFELEKRLTLEPKPDTDKFLETCLEKMQKDTSAGKESILPALLHKPCVPAVSRTDHIWHYDEEYLPDASKPASANTPEQADGGGSNGFIAVNLSSCKQEIDSKEWVIVDKEQDLQDFRTNEAVGHKTTGSPSDEEPEVLQVLEASPQDEKLQLGPWAENDHLKKETSGVVLALSAEGPPTAASEQYTDRLELQPGAASQFIAATPTSLMEAQAEGPLTAITIPRPSVASTQSTSGSFHCGQQPEKKDLQPMEPTVELYSPRENFSGLVVTEGEPPSGGSRTDLGLQIDHIGHDMLPNIRESNKSQDLGPKELPDHNRLVVREFENLPGETEEKSILLESDNEDEKLSRGQHCIEISSLPGDLVIVEKDHSATTEPLDVTKTQTFSVVPNQDKNNEIMKLLTVGTSEISSRDIDPHVEGQIGQVAEMQKNKISKDDDIMSEDLPGHQGDLSTFLHQEGKREKITPRNGELFHCVSENEHGAPTRKDMVRSSFVTRHSRIPVLAQEIDSTLESSSPVSAKEKLLQKKAYQPDLVKLLVEKRQFKSFLGDLSSASDKLLEEKLATVPAPFCEEEVLTPFSRLTVDSHLSRSAEDSFLSPIISQSRKSKIPRPVSWVNTDQVNSSTSSQFFPRPPPGKPPTRPGVEARLRRYKVLGSSNSDSDLFSRLAQILQNGSQKPRSTTQCKSPGSPHNPKTPPKSPVVPRRSPSASPRSSSLPRTSSSSPSRAGRPHHDQRSSSPHLGRSKSPPSHSGSSSSRRSCQQEHCKPSKNGLKGSGSLHHHSASTKTPQGKSKPASKLSR"

nchar(protein)

## TTBK2_Hsap 
##       1244

plist <- list(domain$TTBK2_Hsap, protein$TTBK2_Hsap)
psimmat <- parSeqSim(plist, cores = 4, type = "local", submat = "BLOSUM62")
psimmat

##           [,1]      [,2]
## [1,] 1.0000000 0.4566822
## [2,] 0.4566822 1.0000000

seqalign <- twoSeqSim(domain$TTBK2_Hsap, protein$TTBK2_Hsap)
summary(seqalign)

## [1] "PairwiseAlignmentsSingleSubject object of length 1 with 0 metadata columns"

seqalign

## Local PairwiseAlignmentsSingleSubject (1 of 1)
## pattern:  [1] WKVLRKIGGGGFGEIYDALDMLTRENVALKVE...LMLKHLPPEFSIFLDHISSLDYFTKPDYQLL
## subject: [21] WKVLRKIGGGGFGEIYDALDMLTRENVALKVE...LMLKHLPPEFSIFLDHISSLDYFTKPDYQLL
## score: 1350

#install.packages("doParallel")