Velvet 출력으로부터 N50 계산

N50 길이 정의 및 R로 계산하는 방법 시퀀싱 서클에서 N50 length은 어셈블리 품질을 판단할 때 유용한 휴리스틱입니다. 아래는 N50 length의 정의이며 직관적으로 느껴질 수도, 아닐 수도 있습니다. N50 length은 같은 길이이거나 더 긴 contig들의 합이 전체 contig 길이의 최소 50%가 되는 가장 짧은 contig의 길이입니다. 예시를 위해 어셈블러가 다음과 같은 길이의 contig를 만들었다고 가정해 봅시다 (내림차순): 91 77 70 69 62 56 45 29 16 4 이들의 합은 519bp이므로 N50 length을 만족하는 contig들의 합은 최소 519/2 또는 259.5가 되어야 합니다. Brute force로 계산해 보면 91+77=168, 91+77+70=238, 91+77+70+69=307(충분합니다). 따라서 이 어셈블리의 N50 length은 69bp입니다. 다른 말로는, 이 어셈블리의 핵심 염기서열의 절반 이상이 69bp 이상의 contig에 속한다는 뜻입니다. N50 vs N50 length 기술적으로 N50은 N50 length과는 달리, 최종적으로 한계를 넘어서는 contig의 순서를 의미합니다—이 예에서는 4번째 (69bp가 4번째로 큰 contig)입니다. 불행히도, 더 높은 N50은 더 긴 N50 length와 반대의 의미를 가질 수 있습니다. 일부 논문은 N50 length를 L50이라고 부르며, 대부분은 “length”를 빼고 N50이라고 부르는 편이 간단합니다. 혼란을 최소화하기 위해 N50에 단위를 포함하는 것이 중요합니다. Contig N50 vs Scaffold N50 또 다른 구분은 contig N50과 scaffold N50 사이에 존재합니다. Contig는 연속된 세그먼트이며, scaffold(또는 supercontig)는 gap으로 구분된 contig들로 구성됩니다. Scaffold는 리드 수준에서 paired-end 정보와 대규모 시퀀싱 프로젝트에서는 paired BAC ends를 사용해 구성됩니다. scaffold 시퀀스가 다양한 양의 빈 N으로 채워져 있으므로 scaffold N50은 어셈블리 품질을 비교하는 데에만 사용해서는 안 됩니다. Velvet는 sane expCov 설정과 함께 사용하면 scaffold에 대해 매우 보수적입니다—이로 인해 contigs.fa의 N50은 사실상 contig N50으로 간주될 수 있으며, scaffold N50은 아닙니다. 더 공격적인 프로그램(예: SOAPdenovo)은 별도의 contig 및 scaffold 파일을 생성합니다. Velvet N50 from stats.txt Velvet은 짧은 시퀀스를 위한 인기가 많은 어셈블러이며, Repetitive Align‑Consensus‑Align 접근이 아니라 DeBruijn 그래프와 Eulerian graph 이론을 사용합니다. 어셈블리 과정 중에 N50을 반환하지만, 저는 이 값을 직접 contig에서 유도하고자 했습니다. 이 contig들은 stats.txt라는 표에 요약되어 있습니다. R과 그 누적 합 함수(cumsum)를 사용하면 N50을 쉽게 계산할 수 있습니다. 전략: contig를 크기 내림차순으로 정렬한 뒤, 누적 합이 전체 합의 최소 절반이 되는 첫 번째 값을 찾습니다. 아래는 그 솔루션입니다. 감사합니다: Barry Rowlingson myStatsTablecontigsn50= sum(contigs)/2][1] 주의: kmer에 대한 내용 stats.txt 파일에서 길이는 length=lgth+kmer-1로