2

ここであなたのRスキルが本当に必要です。数日間、このプロットに取り組んできました。私はRの初心者なので、それで説明できるかもしれません。

染色体のシーケンス カバレッジ データがあります (基本的には、すべての染色体の長さに沿った各位置の値であり、ベクトルの長さは数百万になります)。読み取りの素敵なカバレッジ プロットを作成したいと考えています。これは私がこれまでに得たものです： ここに画像の説明を入力

問題ないように見えますが、y ラベルが欠落しているため、どの染色体かがわかります。また、x 軸の変更に問題があったため、カバレッジが終了したところで終了します。さらに、私自身のデータははるかに大きく、特にこのプロットには非常に時間がかかります。これが、この HilbertVis plotLongVector を試した理由です。動作しますが、x 軸、ラベル、y 軸をログに記録する方法、およびベクトルが同じ長さでなくてもプロット上ですべて同じ長さになる方法を理解できません。

source("http://bioconductor.org/biocLite.R")
biocLite("HilbertVis")
library(HilbertVis)
chr1 <- abs(makeRandomTestData(len=1.3e+07)) 
chr2 <- abs(makeRandomTestData(len=1e+07)) 

par(mfcol=c(8, 1), mar=c(1, 1, 1, 1), ylog=T)

# 1st way of trying with some code I found on stackoverflow
# Chr1
plotCoverage <- function(chr1, start, end) { # Defines coverage plotting function.
  plot.new()
  plot.window(c(start, length(chr1)), c(0, 10))
  axis(1, labels=F) 
  axis(4)
  lines(start:end, log(chr1[start:end]), type="l")
}
plotCoverage(chr1, start=1, end=length(chr1)) # Plots coverage result.

# Chr2
plotCoverage <- function(chr2, start, end) { # Defines coverage plotting function.
  plot.new()
  plot.window(c(start, length(chr1)), c(0, 10))
  axis(1, labels=F) 
  axis(4)
  lines(start:end, log(chr2[start:end]), type="l")
}
plotCoverage(chr2, start=1, end=length(chr2)) # Plots coverage result.


# 2nd way of trying with plotLongVector
plotLongVector(chr1, bty="n", ylab="Chr1") # ylab doesn't work
plotLongVector(chr2, bty="n")

それから、特別に関心のある遺伝子と呼ばれる別のベクトルがあります。それらは染色体ベクトルとほぼ同じ長さですが、私のデータでは、値よりも多くのゼロが含まれています。

genes_chr1 <- abs(makeRandomTestData(len=1.3e+07)) 
genes_chr2 <- abs(makeRandomTestData(len=1e+07))

これらの遺伝子ベクターは、染色体の下に赤い点としてプロットしたいと思います! 基本的に、ベクトルに値がある場合 (>0)、長いベクトル プロットの下に点 (または線) として表示されます。これを追加する方法がわかりません！しかし、それはかなり簡単に思えます。

私を助けてください！どうもありがとう。

4

2 に答える 2

4

免責事項:このコードを単純にコピー アンド ペーストして、染色体の位置全体を削除しないでください。位置をサンプリングして (たとえば、@Gx1sptDTDa が示すように)、それらをプロットしてください。そうしないと、コンピューターがドレインに耐えた場合、おそらく何時間も経った後に巨大な黒く塗りつぶされた長方形が表示されます。

を使用するggplot2と、これは非常に簡単に実現できgeom_areaます。ここでは、例を示すために、300 の位置を持つ 3 つの染色体のランダム データを生成しました。あなたはこれに基づいて構築することができます、私は願っています.

# construct a test data with 3 chromosomes and 100 positions
# and random coverage between 0 and 500
set.seed(45)
chr <- rep(paste0("chr", 1:3), each=100)
pos <- rep(1:100, 3)
cov <- sample(0:500, 300)
df  <- data.frame(chr, pos, cov)

require(ggplot2)
p <- ggplot(data = df, aes(x=pos, y=cov)) + geom_area(aes(fill=chr))
p + facet_wrap(~ chr, ncol=1)

ggplot2_geom_area_coverage_plot

于 2013-01-31T17:38:15.513 に答える
1

ggplot2 パッケージを使用できます。

あなたが何を望んでいるのか正確にはわかりませんが、ここで私がしたことは次のとおりです。 ここに画像の説明を入力 これには7000のランダムデータポイントがあります（実際には第1染色体の遺伝子の量の約2倍）. アルファを使用して密集した領域を示しました (ランダム データであるため、ここでは多くはありません)。

library(ggplot2)
Chr1_cov <- sample(1.3e+07,7000)
Chr1 <- data.frame(Cov=Chr1_cov,fil=1)
pl <- qplot(Cov,fil,data=Chr1,geom="pointrange",ymin=0,ymax=1.1,xlab="Chromosome 1",ylab="-",alpha=I(1/50))
print(pl)

以上です。これは 1 秒もかからずに実行されました。ggplot2 には膨大な量の設定があるので、いくつか試してみてください。ファセットを使用して複数のグラフを作成します。

以下のコードは一種の移動平均のためのもので、その出力をプロットしています。実際の移動平均は元の移動平均と (ほぼ) 同じ量のデータ ポイントを持つため、実際の移動平均ではありません。データがより滑らかになるだけです。ただし、このコードは n ポイントごとに平均を取ります。もちろん、かなり高速に実行されますが、多くの詳細情報が失われます。

VeryLongVector <- sample(500,1e+07,replace=TRUE)

movAv <- function(vector,n){
    chops <- as.integer(length(vector)/n)
    count <- 0
    pos <- 0
    Cov <-0
    pos[1:chops] <- 0
    Cov[1:chops] <- 0
    for(c in 1:chops){
        tmpcount <- count + n
        tmppos <- median(count:tmpcount)
        tmpCov <- mean(vector[count:tmpcount])
        pos[c] <- tmppos
        Cov[c] <- tmpCov
        count <- count + n
    }

    result <- data.frame(pos=pos,cov=Cov)
    return(result)
}

Chr1 <- movAv(VeryLongVector,10000)
qplot(pos,cov,data=Chr1,geom="line")

ここに画像の説明を入力

于 2013-01-31T16:38:51.340 に答える