0

異なる染色体上に、数値を分離する要因と見なすことができる数値 (遺伝データ) のデータフレームがあります。次のようになります (各位置のサンプル情報を含む隣接する列があります)。

awk '{print $2 "\t"   $3}' log_values | head 
Chr   Start    Sample1     Sample2
1   102447376  0.46957632  0.38415043
1   102447536  0.49194950  0.30094824
1   102447366  0.49874880 -0.17675325
2   102447366 -0.01910729  0.20264680
1   108332063 -0.03295081  0.07738970
1   109472445  0.02216355 -0.02495788

私がしたいのは、このファイルの他の列から値を取得する一連のプロットを作成することです。行ごとに 1 つをプロットする (別の地域や別のサンプルの結果を表す) 代わりに、[開始] 列の値が互いに十分に近い場合、範囲をカバーするプロットを描画したいと考えています。まず、Start 列に 3 つの値が互いに 1000 以内にある場合にプロットを作成したいと思います。つまり、A から B から C までを含む 1000 であるため、A から B <= 1000 および B から C は <= 1000 ですが、A から C は <= 1000 である必要はありません。以下のコードでは、この 1000 は " CNV_サイズ」。「flanking_size」変数は、プロットを少しズームアウトして、コンテキストを与えるためのものです。

サンプル値を取得すると、Sample1 の 1 つのプロットとして行 1 2 および 3 が強調表示されます。これらのサンプル数は log2Ratios なので、重要なものだけをプロットしたいと思います。私はこれを 0.4 以上または -0.6 以下と定義しています。これは、同じ 3 つの行がサンプル 2 のプロットを生成しないことを意味します。

Chr の列番号/係数が異なるため、4 行目は含まれません。これは、この条件を満たす行の値のみを示す各列の個別のプロットです。したがって、サンプルごとに複数のプロットを作成できますが、この基準を満たす領域の各セットはすべてのサンプルにプロットされます。これが意味をなさない場合は、おそらく、以下の効果のない試みが、私が困惑していることを説明するのに役立つでしょう.

pdf("All_CNVs_closeup.pdf")

CNV_size <- 1000 # bp 
flanking_size <- 1000 # bp

#for(chr in 1:24){
for(chr in 1:1){
 #for(array in 1:24) {
for(array in 1:4) {

 dat <- subset(file, file$Chr == chr ) 
 dat <- subset(dat, dat[,array+6] > 0.4 | dat[,array+6] < -0.6) 
 if(length(dat$Start) > 1 ) { 
 dat <- dat[with(dat, order(Start)), ]

x=dat$Start[2:length(dat$Start)]-dat$Start[1:(length(dat$Start)-1)]
cnv <- 1
while(cnv <= length(x)) { 
for(i in cnv:length(x) ) {
    if(x[i] >= CNV_size) {
        plot_title <- paste(sample_info$Sample.ID[array], files[array],   sep = "   ") 
        plot(dat$Start, -dat[,array+6], main = plot_title , ylim = c(-2,2),      xlim = c(dat$Start[cnv] - flanking_size , dat$Start[i ] + flanking_size) , xlab = chr, ylab = "Log2 Ratio") 
        abline(h = 0.4, col="blue") 
        abline(h = 0, col="red") 
        abline(h = -0.6, col="blue") 
    break
    } # if(x[i] >= CNV_size) {      
    #if(x[i] < CNV_size) i <- i + 1
}   # for(i in cnv:length(x) ) {

cnv <- i 

} # while(x[cnv] <= length(x)) { 
  } # if(length(dat$Start) > 1 ) { 
  } # for(array in 1:24) {
} # for(chr in 1:24){



 dev.off()
4

1 に答える 1