家族のジェノタイプデータを適当に作る

HigashiNopponLikelihood RandomGenotypeFamily
  • 予定パッケージ名
    • HigashiNopponLikelihood
  • 関数名
    • RandomGenotypeFamily
  • タイトル
RandomGenotypeFamily
  • 説明
家系図を与え、指定マーカーについて、アレルとそのアレル頻度情報から、適当にその家系図を満足するジェノタイプデータを作る
  • 使用例
G<-RandomGenotypeFamily(p,NL,Alleles,Probs)
G2<-RandomGenotypeFamily(p,NL,Alleles,Probs,missing=TRUE)
  • ソース
RandomGenotypeFamily<-function(d,n,Alleles,Probs,missing=FALSE){
	ns<-length(d[,1])
	g<-array(0,c(ns,2,n))
	noparents<-which(d[,2]==0,d[,3]==0)
	for(i in noparents){
		for(j in 1:n){
			g[i,,j]<-sample(Alleles[[j]],2,replace=TRUE,Probs[[j]])
		}
	}
	given<-rep(0,ns)
	given[noparents]<-1
	loop<-TRUE
	if(sum(given)==ns)loop<-FALSE
	while(loop){
		yet<-which(given==0)
		for(i in yet){
			prt1<-d[i,2]
			prt2<-d[i,3]
			if(given[prt1]+given[prt2]==2){
				for(j in 1:n){
					g[i,1,j]<-sample(g[prt1,,j],1)
					g[i,2,j]<-sample(g[prt2,,j],1)
				}
				given[i]<-1
			}
		}
		if(sum(given)==ns)loop<-FALSE
	}
	if(missing){
		Type<-p[,5]
		
		G2<-array(0,c(length(p[,1]),2,NL))
		
		G2[which(Type==1),,]<-G[which(Type==1),,]
		g<-G2
	}
	g
}
  • Rdファイル
\name{RandomGenotypeFamily}
\alias{RandomGenotypeFamily}
%- Also NEED an '\alias' for EACH other topic documented here.
\title{
RandomGenotypeFamily
}
\description{
家系図を与え、指定マーカーについて、アレルとそのアレル頻度情報から、適当にその家系図を満足するジェノタイプデータを作る。シミュレーションのため。
}
\usage{
G<-RandomGenotypeFamily(p,NL,Alleles,Probs)
G2<-RandomGenotypeFamily(p,NL,Alleles,Probs,missing=TRUE)
}
\arguments{
\item{p}{家系情報の行列}
\item{Np}{人数}
\item{NL}{マーカー数}
\item{Alleles}{アレル名のベクトルのリスト}
\item{Probs}{アレル頻度のベクトルのリスト}
\item{missing}{TRUEまたはFALSE。TRUEの場合には、DNAサンプル提供者のみジェノタイプがある}
}
%- maybe also 'usage' for other objects documented here.
\details{
%%  ~~ If necessary, more details than the description above ~~
}
\value{
ディプロタイプを表す、３次元アレイ(人数ｘ２ｘマーカー数)
}
\references{
%% ~put references to the literature/web site here ~
}
\author{
%%  ~~who you are~~
}
\note{
%%  ~~further notes~~
}

%% ~Make other sections like Warning with \section{Warning }{....} ~

\seealso{
}

\examples{
p<-matrix(
c(1:8,
  0,0,1,0,3,3,3,3,
  0,0,2,0,4,4,4,4,
  1,0,1,0,1,1,0,0,
  1,1,2,1,1,1,2,1),
  ncol=5)
NL<-15
AandP<-MakeAlleleProb(NL=NL)
Alleles<-AandP$Alleles
Probs<-AandP$Probs
G<-RandomGenotypeFamily(p,NL,Alleles,Probs)
G2<-RandomGenotypeFamily(p,NL,Alleles,Probs,missing=TRUE)
}

家族尤度

HigashiNopponLikelihood RandomGenotype
  • 予定パッケージ名
    • HigashiNopponLikelihood
  • 関数名
    • RandomGenotype
  • タイトル
RandomGenotype
  • 説明
指定マーカーについて、アレルとそのアレル頻度情報から、HWEに基づき、指定人数分ぷジェノタイプを作成する
  • 使用例
Gpool<-RandomGenotype(Np,NL,Alleles,Probs)
  • ソース
RandomGenotype<-function(N,n,Alleles,Probs){
	ns<-N
	g<-array(0,c(ns,2,n))
	noparents<-1:ns
	for(i in noparents){
		for(j in 1:n){
			g[i,,j]<-sample(Alleles[[j]],2,replace=TRUE,Probs[[j]])
		}
	}
	g
}
  • Rdファイル
\name{RandomGenotype}
\alias{RandomGenotype}
%- Also NEED an '\alias' for EACH other topic documented here.
\title{
RandomGenotype
}
\description{
指定マーカーについて、アレルとそのアレル頻度情報から、HWEに基づき、指定人数分ぷジェノタイプを作成する。シミュレーションのため。
}
\usage{
Gpool<-RandomGenotype(Np,NL,Alleles,Probs)
}
\arguments{
\item{Np}{人数}
\item{NL}{マーカー数}
\item{Alleles}{アレル名のベクトルのリスト}
\item{Probs}{アレル頻度のベクトルのリスト}
}
%- maybe also 'usage' for other objects documented here.
\details{
%%  ~~ If necessary, more details than the description above ~~
}
\value{
ディプロタイプを表す、３次元アレイ(人数ｘ２ｘマーカー数)
}
\references{
%% ~put references to the literature/web site here ~
}
\author{
%%  ~~who you are~~
}
\note{
%%  ~~further notes~~
}

%% ~Make other sections like Warning with \section{Warning }{....} ~

\seealso{
}

\examples{
NL<-15
AandP<-MakeAlleleProb(NL=NL)
Alleles<-AandP$Alleles
Probs<-AandP$Probs
Np<-3
Gpool<-RandomGenotype(Np,NL,Alleles,Probs)
}

家族尤度

HigashiNopponLikelihood MakeAlleleProb
  • 予定パッケージ名
    • HigashiNopponLikelihood
  • 関数名
    • MakeAlleleProb
  • タイトル
MakeAlleleProb
  • 説明
指定マーカー数のアレルリストとそのアレル頻度リストを適当に作る
プログラムの検証用のシミュレーションデータ作成のため
  • 使用例
MakeAlleleProb(NL)
  • ソース
MakeAlleleProb<-function(NL=15){
	NL<-15

	Alleles<-NULL
	Probs<-NULL
	library(MCMCpack)
	for(i in 1:NL){
		na<-sample(3:8,1)
		st<-sample(10:20,1)
		tmp<-NULL
		for(j in 1:na){
			a<-paste("",st+j,sep="")
			tmp<-c(tmp,a)
		}
		Alleles[[i]]<-tmp
		Probs[[i]]<-c(rdirichlet(1,rep(1,na)))
		#print(sum(Probs[[i]]))
	}
	list(Alleles=Alleles,Probs=Probs)
}
  • Rdファイル
\name{MakeAlleleProb}
\alias{MakeAlleleProb}
%- Also NEED an '\alias' for EACH other topic documented here.
\title{
MakeAlleleProb
}
\description{
指定マーカー数のアレルリストとそのアレル頻度リストを適当に作る
プログラムの検証用のシミュレーションデータ作成のため
}
\usage{
MakeAlleleProb(NL)
}
\arguments{
自然数
}
%- maybe also 'usage' for other objects documented here.
\details{
%%  ~~ If necessary, more details than the description above ~~
}
\value{
２要素リスト。
第１要素はマーカー数の長さのアレルの名前のリスト。個々のリストは文字列のベクトル。
第２要素はマーカー数の長さのアレルの頻度リスト。個々のリストは数値のベクトルで、その値は0以上かつ、和が１。
}
\references{
%% ~put references to the literature/web site here ~
}
\author{
%%  ~~who you are~~
}
\note{
%%  ~~further notes~~
}

%% ~Make other sections like Warning with \section{Warning }{....} ~

\seealso{
}

\examples{
NL<-15
AandP<-MakeAlleleProb(NL=NL)
Alleles<-AandP$Alleles
Probs<-AandP$Probs
}

家族尤度

HigashiNopponLikelihood MakePedigreeFromFamilyInfo
  • 予定パッケージ名
    • HigashiNopponLikelihood
  • 関数名
    • MakePedigreeFromFamilyInfo
  • タイトル
MakePedigreeFromFamilyInfo
  • 説明
５列の家族情報行列からkinshipパッケージのpedigreeオブジェクトを作る
  • 使用例
MakePedigreeFromFamilyInfo(p)
  • ソース
MakePedigreeFromFamilyInfo<-function(p){
	ns<-length(p[,1])
	affected<-status<-rep(1,ns)
	affected[which(p[,5]==3)]<-0
	status[which(p[,5]==1)]<-0
	ptemp<-pedigree(id=p[,1],dadid=p[,3],momid=p[,2],sex=p[,4],affected=affected,status=status)
	if(sum(ptemp$affected)==0)ptemp$affected<-affected
	ptemp
}
  • Rdファイル
\name{MakePedigreeFromFamilyInfo}
\alias{MakePedigreeFromFamilyInfo}
%- Also NEED an '\alias' for EACH other topic documented here.
\title{
MakePedigreeFromFamilyInfo
}
\description{

}
\usage{
akePedigreeFromFamilyInfo(p)
}
\arguments{
pは５列の行列。第１列はID(数値で1,2,...)。第２列は母親ID。母親不明は0。第３列は父親ID。父親不明は0。第４列は性別。0は男,1は女。第５列は解析用タイプ。1はDNAサンプル提供者。2は探されている家系構成員。3は家系構成員ではあるが、DNAサンプルの提供がなく、かつ、探されていない家系構成員。
}
%- maybe also 'usage' for other objects documented here.
\details{
%%  ~~ If necessary, more details than the description above ~~
}
\value{
pedigreeオブジェクト
}
\references{
%% ~put references to the literature/web site here ~
}
\author{
%%  ~~who you are~~
}
\note{
%%  ~~further notes~~
}

%% ~Make other sections like Warning with \section{Warning }{....} ~

\seealso{
}

\examples{
# needs kinship package
p<-matrix(
c(1:8,
  0,0,1,0,3,3,3,3,
  0,0,2,0,4,4,4,4,
  1,0,1,0,1,1,0,0,
  1,1,2,1,1,1,2,1),
  ncol=5)

ptemp<-MakePedigreeFromFamilyInfo(p)
plot(ptemp)
}

２次元用を多次元錐に拡張

R sphere SpherePowerND 
# df次元球の表面積
SphereSurface<-function(df,r=1){
	2*r^df*pi^(df/2)/gamma(df/2)
}
# df次元球の軸と偏角tであるdf-1次元球の表面積
SphereSurfaceAngle<-function(df,t,r=1){
	d<-r*abs(sin(t))
	ret<-0
	if(df==1){
		ret<-rep(2,length(t))
	}else{
		ret<-SphereSurface(df-1,d)
	}
	ret
}
SphereSurfaceAngle(3,pi/2)
SpherePowerND<-
function (df,Kv, Ks, CheckDists, Tiis, nperm = 1000) 
{
    dfin<-df
    df<-2
    
    Ks[which(Ks == 0)] <- 1e-04
    Kss <- matrix(rep(Kv, length(Ks)), nrow = length(Ks), byrow = TRUE) * 
        Ks
    Ntest <- length(Tiis[, 1])
    tmprss<-seq(from=0,to=1,length.out=nperm)*2*pi
    rss <- cbind(cos(tmprss),sin(tmprss))
    ws<-SphereSurfaceAngle(dfin-1,tmprss)
    ws<-ws/sum(ws)

    PowerOut <- matrix(0, length(Ks), length(CheckDists))
    for (ik in 1:length(Ks)) {
        Pii <- Kss[ik, ]
        AccumProbs <- matrix(0, nperm, length(CheckDists))
        for (i in 1:nperm) {
            rs <- matrix(rss[i, ], nrow = 1)
            tmptheta <- acos(cosxy(t(Pii), rs))
            C <- nrm(Pii) * cos(tmptheta)
            D <- nrm(Pii) * sin(tmptheta)
            As <- t(Tiis %*% t(rs))
            Bs <- t(Tiis %*% Pii)
            current <- which(Bs == max(Bs))[1]
            SerialTopT <- c(current)
            currentR <- 0
            SerialR <- c(currentR)
            SerialStat <- c(Bs[current])
            CheckList <- 1:Ntest
            CheckList <- CheckList[which(CheckList != current)]
            while (length(CheckList) > 0) {
                tmpks <- -(Bs[current] - Bs[CheckList])/(As[current] - 
                  As[CheckList])
                tmpList <- which(tmpks > currentR)
                if (length(tmpList) > 0) {
                  current <- CheckList[which(tmpks == min(tmpks[tmpList]))][1]
                  currentR <- min(tmpks[tmpList])
                  SerialTopT <- c(SerialTopT, current)
                  SerialR <- c(SerialR, currentR)
                  SerialStat <- c(SerialStat, As[current] * currentR + 
                    Bs[current])
                }
                CheckList <- CheckList[tmpList]
                CheckList <- CheckList[which(CheckList != current)]
            }
            intMyX <- intMyX3 <- matrix(0, length(CheckDists), 
                length(SerialTopT))
            AnswerPs <- rep(0, length(CheckDists))
            for (j in 1:length(CheckDists)) {
                tmpSerialStat <- SerialStat - CheckDists[j]
                if (tmpSerialStat[1] >= 0) {
                  AnswerPs[j] <- 1
                }
                if (length(tmpSerialStat) > 1) {
                  for (jj in 1:(length(tmpSerialStat) - 1)) {
                    if (tmpSerialStat[jj] == 0) {
                      intMyX[j, jj] <- 1
                      intMyX3[j, jj] <- sign(As[SerialTopT[jj]])
                    }
                    else {
                      if (tmpSerialStat[jj] * tmpSerialStat[jj + 
                        1] < 0) {
                        intMyX[j, jj] <- 1
                        intMyX3[j, jj] <- sign(As[SerialTopT[jj]])
                      }
                    }
                  }
                }
                if (tmpSerialStat[length(tmpSerialStat)] == 0) {
                  intMyX[j, length(tmpSerialStat)] <- 1
                  intMyX3[j, length(tmpSerialStat)] <- sign(As[SerialTopT[length(tmpSerialStat)]])
                }
                else {
                  if (tmpSerialStat[length(tmpSerialStat)] * 
                    As[SerialTopT[length(tmpSerialStat)]] < 0) {
                    intMyX[j, length(tmpSerialStat)] <- 1
                    intMyX3[j, length(tmpSerialStat)] <- sign(As[SerialTopT[length(tmpSerialStat)]])
                  }
                }
            }
            zeros <- which(apply(abs(intMyX), 1, sum) == 0)
            Crossings <- t(outer(CheckDists, Bs[SerialTopT], 
                FUN = "-"))/As[SerialTopT]
            ChiProbs <- pchisq(Crossings^2, dfin, lower.tail = FALSE)
            AnswerPs <- AnswerPs + apply(t(intMyX3) * ChiProbs, 
                2, sum)
            AnswerPs[zeros] <- AnswerPs[zeros] <- (sign(As[SerialTopT[length(SerialTopT)]]) + 
                1)/2
            AnswerPs[which(AnswerPs < 0)] <- 1 + AnswerPs[which(AnswerPs < 
                0)]
            AnswerPs[which(AnswerPs > 1)] <- AnswerPs[which(AnswerPs > 
                1)] - 1
            AccumProbs[i, ] <- AnswerPs*ws[i]
        }
        #PowOut <- apply(AccumProbs, 2, sum)/nperm
        PowOut <- apply(AccumProbs, 2, sum)
        PowerOut[ik, ] <- PowOut
    }
    list(df = dfin, Kv = Kv, Ks = Ks, CheckDists = CheckDists, 
        Tiis = Tiis, nperm = nperm, p.out = PowerOut)
}


rs<-seq(from=0,to=1,length.out=101)
us<-rs*pi/2
Ks<-5
CheckDists<-5
nperm<-1000
dfs<-c(2,3,4,5,10)
power2markerdfs<-matrix(0,length(dfs),length(us))

for(x in 1:length(dfs)){
	df<-dfs[x]
	ppout<-rep(0,length(us))
	Tiis<-matrix(0,2,2)
	Tiis[1,1]<-1
	Kv<-c(1,0)
	for(i in 1:length(us)){
		Tiis[2,1]<-cos(us[i])
		Tiis[2,2]<-sin(us[i])
		
		tmp<-rbind(Tiis,-Tiis)
		spout<-SpherePowerND(df,Kv,Ks,CheckDists,Tiis,nperm)
		ppout[i]<-c(spout$p.out)
	}

	power2markerdfs[x,]<-ppout

}

matplot(t(power2markerdfs2),type="l")

power2markerdfs2<-matrix(0,length(dfs),length(us))

for(x in 1:length(dfs)){
	df<-dfs[x]
	ppout<-rep(0,length(us))
	Tiis<-matrix(0,3,2)
	Tiis[1,1]<-1
	Kv<-c(1,0)
	for(i in 1:length(us)){
		Tiis[2,1]<-cos(us[i])
		Tiis[2,2]<-sin(us[i])
		Tiis[3,1]<-cos(us[i])
		Tiis[3,2]<--sin(us[i])
		
		tmp<-rbind(Tiis,-Tiis)
		spout<-SpherePowerND(df,Kv,Ks,CheckDists,Tiis,nperm)
		ppout[i]<-c(spout$p.out)
	}

	power2markerdfs2[x,]<-ppout

}

matplot(t(power2markerdfs2),type="l")

マーカーの束ごとにパワーを計算

sphere SpherePower2D 分割表
  • 予定パッケージ名
  • 関数名
    • SpherePower2D
  • タイトル
SpherePower2D
  • 説明
フェノタイプ情報のベクトルと、ジェノタイプ情報の行列があるとす
  • 使用例
SpherePower2D(Kv,Ks,CheckDists,Tiis,nperm)
  • ソース
SpherePower2D<-
function (Kv, Ks, CheckDists, Tiis, nperm = 1000) 
{
    df<-2
    Ks[which(Ks == 0)] <- 1e-04
    Kss <- matrix(rep(Kv, length(Ks)), nrow = length(Ks), byrow = TRUE) * 
        Ks
    Ntest <- length(Tiis[, 1])
    tmprss<-seq(from=0,to=1,length.out=nperm)*2*pi
    rss <- cbind(cos(tmprss),sin(tmprss))
    PowerOut <- matrix(0, length(Ks), length(CheckDists))
    for (ik in 1:length(Ks)) {
        Pii <- Kss[ik, ]
        AccumProbs <- matrix(0, nperm, length(CheckDists))
        for (i in 1:nperm) {
            rs <- matrix(rss[i, ], nrow = 1)
            tmptheta <- acos(cosxy(t(Pii), rs))
            C <- nrm(Pii) * cos(tmptheta)
            D <- nrm(Pii) * sin(tmptheta)
            As <- t(Tiis %*% t(rs))
            Bs <- t(Tiis %*% Pii)
            current <- which(Bs == max(Bs))[1]
            SerialTopT <- c(current)
            currentR <- 0
            SerialR <- c(currentR)
            SerialStat <- c(Bs[current])
            CheckList <- 1:Ntest
            CheckList <- CheckList[which(CheckList != current)]
            while (length(CheckList) > 0) {
                tmpks <- -(Bs[current] - Bs[CheckList])/(As[current] - 
                  As[CheckList])
                tmpList <- which(tmpks > currentR)
                if (length(tmpList) > 0) {
                  current <- CheckList[which(tmpks == min(tmpks[tmpList]))][1]
                  currentR <- min(tmpks[tmpList])
                  SerialTopT <- c(SerialTopT, current)
                  SerialR <- c(SerialR, currentR)
                  SerialStat <- c(SerialStat, As[current] * currentR + 
                    Bs[current])
                }
                CheckList <- CheckList[tmpList]
                CheckList <- CheckList[which(CheckList != current)]
            }
            intMyX <- intMyX3 <- matrix(0, length(CheckDists), 
                length(SerialTopT))
            AnswerPs <- rep(0, length(CheckDists))
            for (j in 1:length(CheckDists)) {
                tmpSerialStat <- SerialStat - CheckDists[j]
                if (tmpSerialStat[1] >= 0) {
                  AnswerPs[j] <- 1
                }
                if (length(tmpSerialStat) > 1) {
                  for (jj in 1:(length(tmpSerialStat) - 1)) {
                    if (tmpSerialStat[jj] == 0) {
                      intMyX[j, jj] <- 1
                      intMyX3[j, jj] <- sign(As[SerialTopT[jj]])
                    }
                    else {
                      if (tmpSerialStat[jj] * tmpSerialStat[jj + 
                        1] < 0) {
                        intMyX[j, jj] <- 1
                        intMyX3[j, jj] <- sign(As[SerialTopT[jj]])
                      }
                    }
                  }
                }
                if (tmpSerialStat[length(tmpSerialStat)] == 0) {
                  intMyX[j, length(tmpSerialStat)] <- 1
                  intMyX3[j, length(tmpSerialStat)] <- sign(As[SerialTopT[length(tmpSerialStat)]])
                }
                else {
                  if (tmpSerialStat[length(tmpSerialStat)] * 
                    As[SerialTopT[length(tmpSerialStat)]] < 0) {
                    intMyX[j, length(tmpSerialStat)] <- 1
                    intMyX3[j, length(tmpSerialStat)] <- sign(As[SerialTopT[length(tmpSerialStat)]])
                  }
                }
            }
            zeros <- which(apply(abs(intMyX), 1, sum) == 0)
            Crossings <- t(outer(CheckDists, Bs[SerialTopT], 
                FUN = "-"))/As[SerialTopT]
            ChiProbs <- pchisq(Crossings^2, df, lower.tail = FALSE)
            AnswerPs <- AnswerPs + apply(t(intMyX3) * ChiProbs, 
                2, sum)
            AnswerPs[zeros] <- AnswerPs[zeros] <- (sign(As[SerialTopT[length(SerialTopT)]]) + 
                1)/2
            AnswerPs[which(AnswerPs < 0)] <- 1 + AnswerPs[which(AnswerPs < 
                0)]
            AnswerPs[which(AnswerPs > 1)] <- AnswerPs[which(AnswerPs > 
                1)] - 1
            AccumProbs[i, ] <- AnswerPs
        }
        PowOut <- apply(AccumProbs, 2, sum)/nperm
        PowerOut[ik, ] <- PowOut
    }
    list(df = df, Kv = Kv, Ks = Ks, CheckDists = CheckDists, 
        Tiis = Tiis, nperm = nperm, p.out = PowerOut)
}
  • Rdファイル
\name{SpherePower2D}
\alias{SpherePower2D}
%- Also NEED an '\alias' for EACH other topic documented here.
\title{
SpherePower2D
}
\description{

}
\usage{
SpherePower2D(Kv,Ks,CheckDists,Tiis,nperm)
}
\arguments{

}
%- maybe also 'usage' for other objects documented here.
\details{
%%  ~~ If necessary, more details than the description above ~~
}
\value{

}
\references{
%% ~put references to the literature/web site here ~
}
\author{
%%  ~~who you are~~
}
\note{
%%  ~~further notes~~
}

%% ~Make other sections like Warning with \section{Warning }{....} ~

\seealso{
}

\examples{
}

マーカーの束ごとにパワーを計算

sphere SpherePowerGene 分割表
  • 予定パッケージ名
  • 関数名
    • SpherePowerGene
  • タイトル
SpherePowerGene
  • 説明
フェノタイプ情報のベクトルと、ジェノタイプ情報の行列があるとする。
フェノタイプベクトルから分割表を作り、
ジェノタイプ行列から、テストベクトルを作る。
遺伝子は、１個以上のマーカーを持ち、
その遺伝子に関する検定は、その遺伝子に属するマーカーに関する複数の検定の最大検定統計量であるものとする。
検定閾値を与え、真の関連をマーカーの方向にあるとみなして、パワーを遺伝子単位で推定する
  • 使用例


    

  • ソース
    • 

SpherePowerGene<-function(p,g,ptype=c(0,1),genes,Ks,CheckDists,nperm){
	O<-MakeOTable(p,ptype)
	Sp<-RegularSpherize(O)
	Ngenes<-length(genes[,1])
	Ntobechecked<-sum(genes[,2])-sum(genes[,1])+Ngenes
	df<-length(p)-1

	powouts<-matrix(0,Ntobechecked,length(CheckDists))
	counter<-1
	for(i in 1:Ngenes){
		W<-MakeWTableS(as.matrix(g[,genes[i,1]:genes[i,2]]))
		Tiis<-matrix(0,length(W[,1]),df)
		for(j in 1:length(W[,1])){
			Tiis[j,]<-TestDf1(O,W[j,],Sp)$Tiist
		}
		Tiis<-rbind(Tiis,-Tiis)
		for(j in 1:(length(Tiis[,1])/2)){
			Kv<-Tiis[j,]
			powouts[counter,]<-SpherePower(df,Kv,Ks,CheckDists,Tiis)$p.out
			counter<-counter+1
		}
	}
	#matplot(t(powouts),type="l")
	list(powers=powouts,p=p,g=g,ptype=ptype,genes=genes,Ks=Ks,CheckDists=CheckDists,nperm=nperm)
}



    

  • Rdファイル
    • 

\name{SpherePowerGene}
\alias{SpherePowerGene}
%- Also NEED an '\alias' for EACH other topic documented here.
\title{
SpherePowerGene
}
\description{
フェノタイプ情報のベクトルと、ジェノタイプ情報の行列があるとする。
フェノタイプベクトルから分割表を作り、
ジェノタイプ行列から、テストベクトルを作る。
遺伝子は、１個以上のマーカーを持ち、
その遺伝子に関する検定は、その遺伝子に属するマーカーに関する複数の検定の最大検定統計量であるものとする。
検定閾値を与え、真の関連をマーカーの方向にあるとみなして、パワーを遺伝子単位で推定する
}
\usage{
SpherePowerGene(p=p,g=g3,genes=LocGenes,Ks=Ks,CheckDists=CheckDists,nperm=nperm)
}
\arguments{
\item{p}{サンプル数の長さのベクトル、フェノタイプ}
\item{g}{サンプル数の行数、マーカー数の列数の行列、ジェノタイプ}
\item{genes}{２列の行列。各行は、遺伝子などの長さのある要素に対応し、第１列の値から第２列の値までのマーカーが、その要素に属するものとする}
\item{Ks}{対立仮説の原点からの距離}
\item{CheckDists}{帰無仮説検定の検定閾値(原点からの距離)}
\item{nperm}{推定に用いる乱方向の数}
}
%- maybe also 'usage' for other objects documented here.
\details{
%%  ~~ If necessary, more details than the description above ~~
}
\value{
引数を格納するとともに、遺伝子上のマーカーのパワーが行列状で返る
}
\references{
%% ~put references to the literature/web site here ~
}
\author{
%%  ~~who you are~~
}
\note{
%%  ~~further notes~~
}

%% ~Make other sections like Warning with \section{Warning }{....} ~

\seealso{
}

\examples{
Ns<-1000
Nm<-1000
k<-c(0,1)
p<-sample(k,Ns,replace=TRUE)
num1<-sample(1:(Ns-5),Nm,replace=TRUE)
num1<-sort(num1)

num2<-rep(0,Nm)
for(i in 1:Nm){
	num2[i]<-sample((num1[i]+1):Ns,1)
}
g<-matrix(0,Ns,Nm)
for(i in 1:Nm){
	g[,i]<-c(rep(0,num1[i]),rep(1,num2[i]-num1[i]),rep(2,Ns-num2[i]))
}

Nm2<-500
g2<-matrix(sample(c(0,1,2),Nm2*Ns,replace=TRUE),nrow=Ns)

g3<-cbind(g,g2)
Nm3<-Nm+Nm2
g3<-g3[,sample(1:Nm3)]

image(g3)

Ngenes<-10

LocGenes<-matrix(sort(sample(1:Nm3,Ngenes*2)),ncol=2,byrow=TRUE)

matplot(LocGenes)
Ks<-3
CheckDists<-0:6
nperm<-1000
powerGenes<-SpherePowerGene(p=p,g=g3,genes=LocGenes,Ks=Ks,CheckDists=CheckDists,nperm=nperm)
matplot(t(powerGenes$powers),type="l")
}