seqinr 점도 - 축 변경

Question

데이터 세트 : seq1 및 seq2 (DNA 시퀀스)가 있습니다. 저는 데이터 플롯을 만들고 두 서열을 비교하고 두 서열이 일치하는 곳에 도트를 배치하기를 원했습니다. seqinr의 점을 사용하여이 작업을 수행 할 수 있었지만 축에 시퀀스를 나열하는 것이 불가능하여 어떤 점이 일치하는지 확인할 수 있습니다. 본질적으로 숫자를 시퀀스 문자로 바꾸고 싶습니다.

어쨌든 이것을 할 수 있습니까? 어쩌면 ggplot2를 통해?

seq1 <- c("G","C","T","A","G","T","C","A","G","A","T","C","T","G","A","C","G","C","T","A") 
seq2 <- c("G","A","T","G","G","T","C","A","C","A","T","C","T","G","C","C","G","C") 

그리고 이것이 내가이 그래프 생성 방법은 다음과 같습니다 :

이 내 순서는

dotPlot(seq1, seq2, main = "Dot plot of 2 different sequences 
\nwsize = 4, wstep = 1, nmatch = 3", wsize = 4, wstep = 1, nmatch = 3) 

Answer 1

어떻게 수정 된 버전 (약 일명 "빠른 앤 더러운 해킹 ")? 첫째, R에 다음 코드를 복사하여 붙여 넣기 :

seq1 <- c("G","C","T","A","G","T","C","A","G","A","T","C","T","G","A","C","G","C","T","A") 
seq2 <- c("G","A","T","G","G","T","C","A","C","A","T","C","T","G","C","C","G","C") 
xy<-dotplot(seq1, seq2, wsize = 4, wstep = 1, nmatch = 3) 

은 다음과 플롯 생성한다 :

dotplot<-function (seq1, seq2, wsize = 1, wstep = 1, nmatch = 1, cols = c("white", 
    "black"), xlab = deparse(substitute(seq1)), ylab = deparse(substitute(seq2)), type=2, 
    ...) { 

    cat("This is a modification of the function dotPlot from package seqinr.\n") 

    require(seqinr) 

    if (nchar(seq1[1]) > 1) 
     stop("seq1 should be provided as a vector of single chars") 
    if (nchar(seq2[1]) > 1) 
     stop("seq2 should be provided as a vector of single chars") 
    if (wsize < 1) 
     stop("non allowed value for wsize") 
    if (wstep < 1) 
     stop("non allowed value for wstep") 
    if (nmatch < 1) 
     stop("non allowed value for nmatch") 
    if (nmatch > wsize) 
     stop("nmatch > wsize is not allowed") 
    mkwin <- function(seq, wsize, wstep) { 
     sapply(seq(from = 1, to = length(seq) - wsize + 1, by = wstep), 
      function(i) c2s(seq[i:(i + wsize - 1)])) 
    } 
    wseq1 <- mkwin(seq1, wsize, wstep) 
    wseq2 <- mkwin(seq2, wsize, wstep) 
    if (nmatch == wsize) { 
     xy <- outer(wseq1, wseq2, "==") 
    } 
    else { 
     "%==%" <- function(x, y) colSums(sapply(x, s2c) == sapply(y, 
      s2c)) >= nmatch 
     xy <- outer(wseq1, wseq2, "%==%") 
    } 

    if(type==1) { 
     image(x = seq(from = 1, to = length(seq1), length = length(wseq1)), 
      y = seq(from = 1, to = length(seq2), length = length(wseq2)), 
      z = xy, col = col, xlab = xlab, ylab = ylab, axes=F, ...) 
     box() 
    } 

    colnames(xy)<-wseq2 
    rownames(xy)<-wseq1 

    xy2<-matrix(nrow=length(seq1), ncol=length(seq2), data=FALSE) 
    rownames(xy2)<-seq1 
    colnames(xy2)<-seq2 
    ind<-which(xy, arr.ind=T) 
    xy2[ind]<-TRUE 
    ind<-data.frame(ind, row.names=NULL)  

    res<-data.frame(row=c(), col=c()) 
    for(i in 1:nrow(ind)) { 
     DF<-data.frame(row=seq(from=ind$row[i], to=ind$row[i]+wsize-1), 
         col=seq(from=ind$col[i], to=ind$col[i]+wsize-1)) 
     res<-rbind(res, DF) 
    } 
    xy2[as.matrix(res)]<-TRUE 

    if(type==2) { 
     image(x = seq(from = 1, to = length(seq1)), y = seq(from = 1, to = length(seq2)), z = xy2, col = cols, xlab = xlab, ylab = ylab, axes=F, ...) 
     box() 
     axis(side=1, at=1:length(seq1), labels=seq1) 
     axis(side=2, at=1:length(seq2), labels=seq2, las=1) 
    } 

    out<-list(type1=xy, type2=xy2) 
    return(out) 
} 

그런 다음 사용자가 제공 한 예를 실행 한마디로

을, 두 시퀀스 사이에서 전체 행렬을 생성하기 위해 함수를 약간 확장했습니다. 출력 객체 xy을 검사하면 원본 (type1) 행렬과 펼쳐진 행 (type2)이 표시됩니다. 매우 긴 서열의 경우이 변경은 효율적이지 않으며 축의 뉴클레오티드/아미노산 라벨이 서로 겹칠 것임은 말할 것도 없습니다. 새 인수 type을 사용하여 type1과 type2 사이의 플롯 유형을 변경할 수 있습니다.

Answer 2

두 시퀀스를 비교하는 대신 TraMineR R 패키지의 seqalign 기능을 고려해 볼 수 있습니다.

library(TraMineR) 
seq1 <- c("G","C","T","A","G","T","C","A","G","A","T","C", 
     "T","G","A","C","G","C","T","A") 
## Filling seq2 with "*" to equalize sequence length 
seq2m <- c("G","A","T","G","G","T","C","A","C","A","T","C", 
     "T","G","C","C","G","C","*","*") 

## defining sequence object interpreting "*"s as mising states 
seq <- seqdef(rbind(seq1,seq2m), missing="*") 

## Setting all substitution costs as 2 
cost <- matrix(rep(2,16),4,4) 
diag(cost) <- 0 
cost 

## Comparing sequences 1 and 2 
sa <- seqalign(seq, 1:2, indel=1, sm=cost) 
print(sa) 
plot(sa) 

줄거리 정합 소자 (EQU) 다른으로 하나의 시퀀스를 변환 최소한의 수정 작업 (SUB 및 IND)을 나타낸다 : 여기

는 일례이다.