R의 xtabs를 사용한 총 가중 읽기 커버리지

R의 xtabs를 활용한 총 가중 읽기 커버리지 계산Samtools와 BioPerl 래퍼 Bio::DB::Sam은 일반적으로 각 염기쌍마다 깊이(Depth)를 제공합니다. 이는 최소 한 개의 read가 정렬된 모든 위치에 대해 깊이 값을 담은 배열입니다. 이 방식은 제가 사용하는데는 충분하지만, 트랜스크립트와 같은 타깃 중에서 가장 높은 총 가중 읽기 커버리지를 가진 타깃을 빠르게 찾고 싶습니다. 가중 읽기 커버리지란 각 read가 정렬된 모든 염기쌍의 깊이 합계를 말합니다.이 문제를 해결하기 위해, 타깃-포지션-깊이 정보를 포함하는 테이블을 R에 불러옵니다. 필요한 최소 컬럼은 다음과 같습니다.target (타깃 이름)pos (염기쌍 위치)depth (해당 위치에서의 깊이)예시 테이블(포지션과 깊이만 표시):target pos depth tx500090 227 12 tx500090 228 13 tx500090 229 14 tx500090 230 15 tx500090 231 1 ... tx500123 184 167 tx500123 185 168 tx500123 186 169 tx500123 187 270 tx500123 188 271 tx500123 189 2R에서 xtabs()함수를 사용하면 기본적으로 깊이 값을 가중치로 삼아 포지션을 합산할 수 있습니다. 다음과 같이 사용합니다.myCoverage <- myxTab <- xtabs(~ target + pos + depth, data = myData) head(myxTab)위 코드는 기본적으로 depth를 가중치로 사용해 각 타깃의 총 가중 읽기 커버리지를 계산하고, 정렬되지 않은 결과를 반환합니다. 예시 결과는 다음과 같습니다.target tx500090 tx500123 tx500134 tx500155 tx500170 x 38 92 610 46 176 ...이 방법은 다음과 같은 수식 대신 약 10,000배 빠릅니다.by(myCoverage, myCoverage$tx, function(x) { sum(x$depth) })따라서, xtabs()를 활용하면 가중 읽기 커버리지를 효율적으로 계산할 수 있습니다.