[その他、ソフト]VAST 2010 Mini3の解析

塩基配列の変異と症状との関連を探る。

まずNativeの10種類配列をチェックすると、1404塩基のうち1234塩基は同じで120塩基が異なっていた。

異なっていたというのは、10種のうちどれかが違っているということ。

1234塩基は10種で同じなので、この塩基のどこかに変異があれば怪しいということになる。

次に58患者の配列と比較すると、1234塩基のうち異なっている塩基の数は次の通り。

f:id:isseing333:20100513173454j:image

変異が起きているのは12～16塩基と少ない（順に1,7,9,19,22人）。

次に10Nativeグループで共通な塩基のうち、58患者の中でも違った変異パターンをしている57塩基でヒートマップを描く。

f:id:isseing333:20100513173629j:image

患者が3グループに分かれてるっぽかったので3つに分割し、Symptomと比較してみた。

1 2 3
Mild 4 4 9
Mod 10 4 13
Sev 10 0 4

Severe14人のうち10人がグループ1に居る。

ヒートマップでグループ1（右側のグループ）をみると、842番の塩基があやしい。

842番の塩基がT→Cへと変異している患者は58人中17人居て、そのうち9人がSevereであった（53%）。

842番が変異していない患者のSevere割合は12%であったので、リスク比は4.4倍！！

ちなみに、10Nativeグループで共通で、かつ58患者で同じ変異をしているのは12塩基で、次のような変異の仕方だった。

  • 442 :C→G
  • 677 :T→G
  • 689 :G→C
  • 781 :G→A
  • 895 :C→G
  • 910 :G→A
  • 962 :G→T
  • 1115:A→T
  • 1153:C→T
  • 1328:T→G
  • 1334:T→C
  • 1395:G→T

コードは下記です↓

#---遺伝子配列のデータはcsvに加工しています（1列目が対象者、2列目が塩基配列）
#---------Nativeの遺伝子配列を確認してみる
Native <- read.csv("D:/R/VAST 2010/VAST2010-Mini-3-Data Files-04-09/NativeSequences.csv", header=F)

#------塩基配列をバラす(1塩基目が空欄なので1404塩基ある)
NativeGene <- matrix(rep("a", 1405*nrow(Native)), 1405, nrow(Native))
colnames(NativeGene) <- Native[, 1]
for(i in 1:1405){
	for(j in 1:nrow(Native)){
		NativeGene[i, j] <- substr(as.character(Native[j, 2]), i, i)
	}
}

#------全て同じ塩基のフラグ
SameGene <- ifelse(
	NativeGene[, 1]==NativeGene[, 2] & NativeGene[, 1]==NativeGene[, 3] & 
	NativeGene[, 1]==NativeGene[, 4] & NativeGene[, 1]==NativeGene[, 5] & 
	NativeGene[, 1]==NativeGene[, 6] & NativeGene[, 1]==NativeGene[, 7] & 
	NativeGene[, 1]==NativeGene[, 8] & NativeGene[, 1]==NativeGene[, 9] & 
	NativeGene[, 1]==NativeGene[, 10], 1, 0)
table(SameGene)		#170個がどれかのグループで違う塩基→他の1234は10グループで同じなのだから、そこの塩基が違っていれば問題であろう


#---------患者のデータ
Patient <- read.csv("D:/R/VAST 2010/VAST2010-Mini-3-Data Files-04-09/CurrentOutbreakSequences.csv", header=F)

#------塩基配列をバラす
PatientGene <- matrix(rep("a", 1405*nrow(Patient)), 1405, nrow(Patient))
colnames(PatientGene) <- Patient[, 1]
for(i in 1:1405){
	for(j in 1:nrow(Patient)){
		PatientGene[i, j] <- substr(as.character(Patient[j, 2]), i, i)
	}
}

#------Nativeと結合
PNGene <- data.frame(PatientGene, NativeGene)
SamePN <- matrix(rep(-1, 1405, 58), 1405, 58)
for(i in 1:58) SamePN[, i] <- ifelse(PNGene[, i]==PNGene[, 59], 1, 0)

#------1234塩基のうち異なっている塩基の数（同じ塩基が空欄の分1つ多い、異なっている塩基の数には影響なし）
Ndiff <- apply(SamePN, 2, function(d) table(d[SameGene==1])[1])

library(ggplot2)
qplot(Ndiff)
table(Ndiff)		#12～16個（1,7,9,19,22）

#------58サンプルが同じ塩基かどうか
PatientGeneN <- PatientGene
PatientGeneN[PatientGeneN=="A"] <- 1
PatientGeneN[PatientGeneN=="T"] <- 2
PatientGeneN[PatientGeneN=="G"] <- 3
PatientGeneN[PatientGeneN=="C"] <- 4
SamePP <- apply(PatientGeneN, 1, function(d) ifelse(mean(as.numeric(d))==as.numeric(d[1]), 1, 0))
table(SamePP)		#57個違う
table(SamePP, SameGene)	#このうち56個がNative10グループでは同じ


#------ヒートマップ
library(gplots)

rownames(SamePN) <- 0:1404	#1～1404が塩基の順番を示している
colnames(SamePN) <- Patient[,1]
jpeg("D:/R/VAST 2010/VAST2010-Mini-3-Data Files-04-09/患者で変異のある塩基.jpeg", width=600, height=600)
	heatmap.2(SamePN[SameGene==1 & SamePP==0 & !is.na(SamePP),], col=redgreen(64), 
		density.info="none", trace="none", margins=c(2, 3), 
		hclustfun=function(d) hclust(d, method="ward"))
dev.off()

#------患者が3つのグループに分かれている
plot(PatientTree <- hclust(dist(t(SamePN[SameGene==1 & SamePP==0 & !is.na(SamePP),])), method="ward"))
rect.hclust(PatientTree, k=3)
PatientGroup <- cutree(PatientTree, k=3)

PatientGroup <- data.frame(Sequence_ID=Patient[,1], class=PatientGroup)

#---------患者の背景との関連をチェック
Character <- read.table("D:/R/VAST 2010/VAST2010-Mini-3-Data Files-04-09/DiseaseCharacteristics.txt", header=T)

CharacterGroup <- merge(Character, PatientGroup, by="Sequence_ID")
(Table1 <- table(CharacterGroup$Symptoms, CharacterGroup$class))
chisq.test(Table1)
Table1
round(prop.table(Table1, 1), digits=2)
round(prop.table(Table1, 2), digits=2)


#---------842番があやしい
P842 <- data.frame(name = colnames(SamePN[, SamePN[843, ]==0]), P842 = 1)	#842番に変異がある対象者17人
CharP842 <- merge(Character, P842, by.x="Sequence_ID", by.y="name", all=T)		#17人の背景
CharP842[is.na(CharP842)] <- 0
Table2 <- table(CharP842$Symptom, CharP842$P842)

Table2
round(prop.table(Table2, 1), digits=2)
round(prop.table(Table2, 2), digits=2)


#どのように変異しているのか
PatientGene[843, SamePN[843, ]==0]
NativeGene[843, ]			#T→Cへと変異

ページTOPへ