SciPy 스택으로 게놈에 대한 포괄적인 소개

생물정보학은 생물학적 데이터를 분석하고 이해하기 위한 방법과 소프트웨어 도구를 개발하는 학제 간 분야입니다.

간단히 말해서 생물학을 위한 데이터 과학이라고 생각하시면 됩니다.

많은 종류의 생물학적 데이터 중 게놈 데이터는 가장 널리 분석되는 데이터 중 하나입니다. 특히, 차세대 DNA 염기서열분석(NGS) 기술의 급속한 발전으로 게놈 데이터의 양이 기하급수적으로 증가하고 있습니다. Stephens, Zachary D et al.에 따르면 게놈 데이터 수집은 연간 엑사바이트 규모입니다.

이 포스트에서는 SciPy Stack을 사용하여 인간 게놈에 대한 GFF3 파일을 분석하는 예를 시연합니다. GFF3(Generic Feature Format Version 3)은 게놈 기능을 저장하기 위한 현재 표준 텍스트 파일 형식입니다. 특히, 이 게시물에서는 SciPy 스택을 사용하여 인간 게놈에 대한 다음 질문에 답하는 방법을 배웁니다.

1. 얼마나 많은 게놈이 불완전합니까?
2. 게놈에는 몇 개의 유전자가 있습니까?
3. 전형적인 유전자의 길이는 얼마입니까?
4. 염색체 사이의 유전자 분포는 어떻게 생겼습니까?

인간 게놈에 대한 최신 GFF3 파일은 여기에서 다운로드할 수 있습니다. 이 디렉토리에 있는 README 파일은 이 데이터 형식에 대한 간략한 설명을 제공하며 더 자세한 사양은 여기에서 찾을 수 있습니다.

빠르고 유연하며 표현력이 풍부한 데이터 구조를 제공하는 SciPy 스택의 주요 구성 요소인 Pandas를 사용하여 GFF3 파일을 조작하고 이해합니다.

설정

먼저 SciPy 스택이 설치된 가상 환경을 설정해 보겠습니다. 이 프로세스는 소스에서 수동으로 빌드하는 경우 시간이 많이 소요될 수 있습니다. 스택에는 외부 FORTRAN 또는 C 코드에 따라 달라지는 많은 패키지가 포함되기 때문입니다. 여기서는 설정이 매우 쉬운 Miniconda를 사용하는 것이 좋습니다.

 wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh bash Miniconda3-latest-Linux-x86_64.sh -b

bash 라인의 -b 플래그는 배치 모드에서 실행하도록 지시합니다. 위의 명령을 사용하여 Miniconda를 성공적으로 설치한 후 게놈을 위한 새 가상 환경을 시작한 다음 SciPy 스택을 설치합니다.

 mkdir -p genomics cd genomics conda create -p venv ipython matplotlib pandas

이 게시물에서 사용할 패키지는 3개만 지정했습니다. 모든 패키지를 SciPy 스택에 나열하려면 conda create 명령 끝에 패키지를 추가하기만 하면 됩니다.

패키지의 정확한 이름이 확실하지 않은 경우 conda search 를 시도하십시오. 가상 환경을 활성화하고 IPython을 시작합시다.

 source activate venv/ ipython

IPython은 기본 Python 인터프리터 인터페이스를 훨씬 더 강력하게 대체하므로 기본 Python 인터프리터에서 수행했던 모든 작업을 IPython에서도 수행할 수 있습니다. 아직 IPython을 사용해 본 적이 없는 모든 Python 프로그래머가 시도해 볼 것을 적극 권장합니다.

주석 파일 다운로드

이제 설정이 완료되었으므로 GFF3 형식의 인간 게놈 주석 파일을 다운로드하겠습니다.

약 37MB로, 일반 텍스트로 약 3GB인 인간 게놈의 정보 내용에 비하면 매우 작은 파일이다. GFF3 파일에는 시퀀스의 주석만 포함되어 있는 반면 시퀀스 데이터는 일반적으로 FASTA라는 다른 파일 형식으로 저장되기 때문입니다. 관심이 있는 경우 여기에서 FASTA를 다운로드할 수 있지만 이 자습서에서는 시퀀스 데이터를 사용하지 않습니다.

 !wget ftp://ftp.ensembl.org/pub/release-85/gff3/homo_sapiens/Homo_sapiens.GRCh38.85.gff3.gz

접두어 ! IPython은 이것이 Python 명령 대신 쉘 명령임을 알려줍니다. 그러나 IPython은 접두사 없이도 ls , pwd , rm , mkdir , rmdir 과 같이 자주 사용되는 일부 쉘 명령을 처리할 수 있습니다 ! .

GFF 파일의 헤드를 살펴보면 ## 또는 #! .

README에 따르면 ## 은 메타데이터가 안정적인 반면 #! 실험적이라는 뜻이다.

나중에 사양에 따라 더 미묘한 의미를 가진 또 다른 지시문인 ### 도 볼 수 있습니다.

사람이 읽을 수 있는 주석은 # 하나 뒤에 와야 합니다. 단순화를 위해 # 으로 시작하는 모든 줄을 주석으로 처리하고 분석 중에 무시합니다.

 ##gff-version 3 ##sequence-region 1 1 248956422 ##sequence-region 10 1 133797422 ##sequence-region 11 1 135086622 ##sequence-region 12 1 133275309 ... ##sequence-region MT 1 16569 ##sequence-region X 1 156040895 ##sequence-region Y 2781480 56887902 #!genome-build GRCh38.p7 #!genome-version GRCh38 #!genome-date 2013-12 #!genome-build-accession NCBI:GCA_000001405.22 #!genebuild-last-updated 2016-06

첫 번째 줄은 이 파일에 사용된 GFF 형식의 버전이 3임을 나타냅니다.

그 다음은 모든 시퀀스 영역의 요약입니다. 나중에 살펴보겠지만 이러한 정보는 파일의 본문 부분에서도 찾을 수 있습니다.

#! 이 주석 파일이 적용되는 게놈, GRCh38.p7의 특정 빌드에 대한 정보를 보여줍니다.

GCR(게놈 참조 컨소시엄)은 인간, 마우스, 제브라피시 및 닭을 포함하여 여러 참조 게놈 어셈블리에 대한 업데이트 및 개선을 감독하는 국제 컨소시엄입니다.

이 파일을 살펴보면 처음 몇 줄의 주석이 있습니다.

 1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_000001.11 ### 1 . biological_region 10469 11240 1.3e+03 . . external_name=oe %3D 0.79;logic_name=cpg 1 . biological_region 10650 10657 0.999 + . logic_name=eponine 1 . biological_region 10655 10657 0.999 - . logic_name=eponine 1 . biological_region 10678 10687 0.999 + . logic_name=eponine 1 . biological_region 10681 10688 0.999 - . logic_name=eponine ...

열은 seqid , 소스 , 유형 , 시작 , 종료 , 점수 , 가닥 , 단계 , 속성 입니다. 그들 중 일부는 매우 이해하기 쉽습니다. 첫 번째 줄을 예로 들어 보겠습니다.

 1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_000001.11

이것은 첫 번째 염기부터 24,895,622번째 염기까지 seqid가 1인 첫 번째 염색체의 주석입니다.

즉, 첫 번째 염색체의 길이는 약 2,500만 염기입니다.

우리의 분석에는 값이 인 세 열의 정보가 필요하지 않습니다 . (즉, 점수, 가닥 및 위상), 따라서 지금은 단순히 무시할 수 있습니다.

마지막 속성 열에는 염색체 1에도 CM000663.2, chr1 및 NC_000001.11이라는 세 개의 별칭이 있습니다. 이것이 기본적으로 GFF3 파일의 모양이지만, 한 줄씩 검사하지 않을 것이므로 전체 파일을 Pandas에 로드할 시간입니다.

Pandas는 탭으로 구분된 파일이기 때문에 GFF3 형식을 처리하는 데 적합하고 Pandas는 CSV와 같은 파일 읽기를 매우 잘 지원합니다.

탭으로 구분된 형식의 한 가지 예외는 GFF3에 ##FASTA 가 포함된 경우입니다.

사양에 따르면 ##FASTA 는 주석 부분의 끝을 나타내며 FASTA(탭으로 구분되지 않음) 형식의 하나 이상의 시퀀스가 ​​뒤따릅니다. 그러나 이것은 우리가 분석할 GFF3 파일의 경우가 아닙니다.

 In [1]: import pandas as pd In [2]: pd.__version__ Out[2]: '0.18.1' In [3]: col_names = ['seqid', 'source', 'type', 'start', 'end', 'score', 'strand', 'phase', 'attributes'] Out[3]: df = pd.read_csv('Homo_sapiens.GRCh38.85.gff3.gz', compression='gzip', sep='\t', comment='#', low_memory=False, header=None, names=col_names)

위의 마지막 줄은 pandas.read_csv 메서드를 사용하여 전체 GFF3 파일을 로드합니다.

표준 CSV 파일이 아니므로 호출을 약간 사용자 정의해야 합니다.

먼저, header=None 을 사용하여 GFF3에서 헤더 정보를 사용할 수 없음을 Pandas에 알린 다음 names=col_names 로 각 열의 정확한 이름을 지정합니다.

names 인수를 지정하지 않으면 Pandas는 증분 숫자를 각 열의 이름으로 사용합니다.

sep='\t' 는 Pandas에 열이 쉼표로 구분되는 대신 탭으로 구분된다는 것을 알려줍니다. sep= 에 대한 값은 실제로 정규식(regex)일 수 있습니다. 현재 파일이 각 열에 대해 서로 다른 구분 기호를 사용하는 경우에 편리합니다(예, 그런 일이 발생합니다). comment='#' 는 # 으로 시작하는 행이 주석으로 간주되어 무시됨을 의미합니다.

compression='gzip' 은 Pandas에 입력 파일이 gzip으로 압축되었음을 알려줍니다.

또한 pandas.read_csv 에는 다양한 종류의 CSV와 같은 파일 형식을 읽을 수 있는 풍부한 매개변수 세트가 있습니다.

반환되는 값의 유형은 Pandas에서 가장 중요한 데이터 구조인 DataFrame 으로 2D 데이터를 표현하는 데 사용됩니다.

Pandas는 또한 각각 1D 및 3D 데이터에 대한 Series 및 Panel 데이터 구조를 가지고 있습니다. Pandas의 데이터 구조에 대한 소개는 문서를 참조하십시오.

.head 메소드를 사용하여 처음 몇 개 항목을 살펴보겠습니다.

 In [18]: df.head() Out[18]: seqid source type start end score strand phase attributes 0 1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_00000... 1 1 . biological_region 10469 11240 1.3e+03 . . external_name=oe %3D 0.79;logic_name=cpg 2 1 . biological_region 10650 10657 0.999 + . logic_name=eponine 3 1 . biological_region 10655 10657 0.999 - . logic_name=eponine 4 1 . biological_region 10678 10687 0.999 + . logic_name=eponine

출력은 속성 열의 더 긴 문자열이 부분적으로 ... 로 대체된 표 형식으로 멋지게 지정됩니다.

pd.set_option('display.max_colwidth', -1) 을 사용하여 긴 문자열을 생략하지 않도록 Pandas를 설정할 수 있습니다. 또한 Pandas에는 사용자 정의할 수 있는 많은 옵션이 있습니다.

다음으로 .info 메서드를 사용하여 이 데이터 프레임에 대한 몇 가지 기본 정보를 알아보겠습니다.

 In [20]: df.info() <class 'pandas.core.frame.DataFrame'> RangeIndex: 2601849 entries, 0 to 2601848 Data columns (total 9 columns): seqid object source object type object start int64 end int64 score object strand object phase object attributes object dtypes: int64(2), object(7) memory usage: 178.7+ MB

이것은 GFF3에 2,601,848개의 주석이 달린 행이 있고 각 행에 9개의 열이 있음을 보여줍니다.

각 열에 대해 해당 데이터 유형도 표시됩니다.

시작과 끝은 int64 유형이며, 게놈의 위치를 ​​나타내는 정수입니다.

다른 열은 모두 object 유형이며, 이는 해당 값이 정수, 부동 소수점 및 문자열의 혼합으로 구성되어 있음을 의미합니다.

모든 정보의 크기는 메모리에 저장된 약 178.7MB 이상입니다. 이것은 약 402MB가 될 압축되지 않은 파일보다 더 작은 것으로 판명되었습니다. 빠른 확인은 아래와 같습니다.

 gunzip -c Homo_sapiens.GRCh38.85.gff3.gz > /tmp/tmp.gff3 && du -s /tmp/tmp.gff3 402M /tmp/tmp.gff3

높은 수준의 보기에서 전체 GFF3 파일을 Python의 DataFrame 개체에 로드했으며 다음 분석은 모두 이 단일 개체를 기반으로 합니다.

이제 첫 번째 열 seqid가 무엇인지 살펴보겠습니다.

 In [29]: df.seqid.unique() Out[29]: array(['1', '10', '11', '12', '13', '14', '15', '16', '17', '18', '19', '2', '20', '21', '22', '3', '4', '5', '6', '7', '8', '9', 'GL000008.2', 'GL000009.2', 'GL000194.1', 'GL000195.1', ... 'KI270757.1', 'MT', 'X', 'Y'], dtype=object) In [30]: df.seqid.unique().shape Out[30]: (194,)

df.seqid 는 데이터 프레임에서 열 데이터에 액세스하는 한 가지 방법입니다. 또 다른 방법은 더 일반적인 구문인 df['seqid'] 입니다. 왜냐하면 열 이름이 Python 예약어(예: class )이거나 . 또는 공백 문자가 있으면 첫 번째 방법( df.seqid )이 작동하지 않습니다.

출력은 염색체 1~22, X, Y 및 미토콘드리아(MT) DNA를 포함하는 194개의 고유한 seqid와 169개의 다른 seqid가 있음을 보여줍니다.

KI 및 GL로 시작하는 seqid는 염색체에 성공적으로 조립되지 않은 게놈의 DNA 서열 또는 스캐폴드입니다.

유전체학에 익숙하지 않은 사람들에게는 이것이 중요합니다.

최초의 인간 게놈 초안이 나온 지 15년이 넘었지만 현재 인간 게놈은 여전히 ​​불완전합니다. 이러한 서열을 조립하는 데 어려움이 있는 것은 주로 게놈의 복잡한 반복 영역 때문입니다.

다음으로 소스 열을 살펴보겠습니다.

README는 소스가 이 기능을 생성한 알고리즘 또는 운영 절차를 설명하기 위한 자유 텍스트 한정자라고 말합니다.

 In [66]: df.source.value_counts() Out[66]: havana 1441093 ensembl_havana 745065 ensembl 228212 . 182510 mirbase 4701 GRCh38 194 insdc 74

이것은 범주형 변수의 빠른 계산에 매우 유용한 value_counts 메서드 사용의 예입니다.

결과에서 이 열에 7개의 가능한 값이 있으며 GFF3 파일의 항목 대부분은 havana, ensembl 및 ensembl_havana에서 가져온 것임을 알 수 있습니다.

이 게시물에서 이러한 소스의 의미와 소스 간의 관계에 대해 자세히 알아볼 수 있습니다.

일을 단순하게 유지하기 위해 GRCh38, havana, ensembl 및 ensembl_havan.a 소스의 항목에 중점을 둘 것입니다.

얼마나 많은 게놈이 불완전합니까?

각 전체 염색체에 대한 정보는 소스 GRCh38의 항목에 있으므로 먼저 나머지를 필터링하고 필터링된 결과를 새 변수 gdf 에 할당하겠습니다.

 In [70]: gdf = df[df.source == 'GRCh38'] In [87]: gdf.shape Out[87]: (194, 9) In [84]: gdf.sample(10) Out[84]: seqid source type start end score strand phase attributes 2511585 KI270708.1 GRCh38 supercontig 1 127682 . . . ID=supercontig:KI270708.1;Alias=chr1_KI270708v... 2510840 GL000208.1 GRCh38 supercontig 1 92689 . . . ID=supercontig:GL000208.1;Alias=chr5_GL000208v... 990810 17 GRCh38 chromosome 1 83257441 . . . ID=chromosome:17;Alias=CM000679.2,chr17,NC_000... 2511481 KI270373.1 GRCh38 supercontig 1 1451 . . . ID=supercontig:KI270373.1;Alias=chrUn_KI270373... 2511490 KI270384.1 GRCh38 supercontig 1 1658 . . . ID=supercontig:KI270384.1;Alias=chrUn_KI270384... 2080148 6 GRCh38 chromosome 1 170805979 . . . ID=chromosome:6;Alias=CM000668.2,chr6,NC_00000... 2511504 KI270412.1 GRCh38 supercontig 1 1179 . . . ID=supercontig:KI270412.1;Alias=chrUn_KI270412... 1201561 19 GRCh38 chromosome 1 58617616 . . . ID=chromosome:19;Alias=CM000681.2,chr19,NC_000... 2511474 KI270340.1 GRCh38 supercontig 1 1428 . . . ID=supercontig:KI270340.1;Alias=chrUn_KI270340... 2594560 Y GRCh38 chromosome 2781480 56887902 . . . ID=chromosome:Y;Alias=CM000686.2,chrY,NC_00002...

Pandas에서는 필터링이 쉽습니다.

df.source == 'GRCh38' 식에서 평가된 값을 검사하면 df 와 동일한 인덱스를 가진 각 항목에 대한 일련의 True 및 False 값입니다. df[] 에 전달하면 해당 값이 True인 항목만 반환됩니다.

df.source == 'GRCh38' 인 df[] 에는 194개의 키가 있습니다.

이전에 보았듯이 seqid 열에는 194개의 고유 값도 있습니다. 이는 gdf의 각 항목이 특정 gdf 에 해당함을 의미합니다.

그런 다음 sample 방법으로 10개의 항목을 무작위로 선택하여 자세히 살펴봅니다.

조립되지 않은 시퀀스는 supercontig 유형이고 다른 시퀀스는 염색체임을 알 수 있습니다. 불완전한 게놈의 부분을 계산하려면 먼저 전체 게놈의 길이를 알아야 합니다. 이 길이는 모든 seqid의 길이를 합한 것입니다.

 In [90]: gdf = gdf.copy() In [91]: gdf['length'] = gdf.end - gdf.start + 1 In [93]: gdf.head() Out[93]: seqid source type start end score strand phase attributes length 0 1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_00000... 248956421 235068 10 GRCh38 chromosome 1 133797422 . . . ID=chromosome:10;Alias=CM000672.2,chr10,NC_000... 133797421 328938 11 GRCh38 chromosome 1 135086622 . . . ID=chromosome:11;Alias=CM000673.2,chr11,NC_000... 135086621 483370 12 GRCh38 chromosome 1 133275309 . . . ID=chromosome:12;Alias=CM000674.2,chr12,NC_000... 133275308 634486 13 GRCh38 chromosome 1 114364328 . . . ID=chromosome:13;Alias=CM000675.2,chr13,NC_000... 114364327 In [97]: gdf.length.sum() Out[97]: 3096629532 In [99]: chrs = [str(_) for _ in range(1, 23)] + ['X', 'Y', 'MT'] In [101]: gdf[-gdf.seqid.isin(chrs)].length.sum() / gdf.length.sum() Out[101]: 0.0037021917421198327

위의 스니펫에서 먼저 .copy() 를 사용하여 gdf 의 복사본을 만들었습니다. 그렇지 않으면 원래 gdf 는 df 의 조각일 뿐이며 직접 수정하면 SettingWithCopyWarning 이 발생합니다(자세한 내용은 여기 참조).

그런 다음 각 항목의 길이를 계산하고 "length"라는 새 열로 gdf 에 다시 추가합니다. 전체 길이는 약 31억으로 밝혀졌고, 조립되지 않은 시퀀스의 비율은 약 0.37%입니다.

다음은 마지막 두 명령에서 슬라이싱이 작동하는 방식입니다.

먼저, 우리는 모두 염색체와 미토콘드리아인 잘 조립된 시퀀스의 모든 seqid를 포함하는 문자열 목록을 만듭니다. 그런 다음 isin 메소드를 사용하여 seqid가 chrs 목록에 있는 모든 항목을 필터링합니다.

목록에 없는 모든 것을 실제로 원하기 때문에(즉, KI 및 GL로 시작하는 어셈블되지 않은 것을 원함) 선택을 반전시키기 위해 인덱스 시작 부분에 빼기 기호( - )가 추가됩니다.

참고: 어셈블된 시퀀스와 어셈블되지 않은 시퀀스는 유형 열로 구분되므로 마지막 줄을 다음과 같이 다시 작성하여 동일한 결과를 얻을 수도 있습니다.

 gdf[(gdf['type'] == 'supercontig')].length.sum() / gdf.length.sum()

얼마나 많은 유전자가 있습니까?

여기서는 대부분의 주석 항목이 속하는 소스 ensembl, havana 및 ensembl_havana의 항목에 중점을 둡니다.

 In [109]: edf = df[df.source.isin(['ensembl', 'havana', 'ensembl_havana'])] In [111]: edf.sample(10) Out[111]: seqid source type start end score strand phase attributes 915996 16 havana CDS 27463541 27463592 . - 2 ID=CDS:ENSP00000457449;Parent=transcript:ENST0... 2531429 X havana exon 41196251 41196359 . + . Parent=transcript:ENST00000462850;Name=ENSE000... 1221944 19 ensembl_havana CDS 5641740 5641946 . + 0 ID=CDS:ENSP00000467423;Parent=transcript:ENST0... 243070 10 havana exon 13116267 13116340 . + . Parent=transcript:ENST00000378764;Name=ENSE000... 2413583 8 ensembl_havana exon 144359184 144359423 . + . Parent=transcript:ENST00000530047;Name=ENSE000... 2160496 6 havana exon 111322569 111322678 . - . Parent=transcript:ENST00000434009;Name=ENSE000... 839952 15 havana exon 76227713 76227897 . - . Parent=transcript:ENST00000565910;Name=ENSE000... 957782 16 ensembl_havana exon 67541653 67541782 . + . Parent=transcript:ENST00000379312;Name=ENSE000... 1632979 21 ensembl_havana exon 37840658 37840709 . - . Parent=transcript:ENST00000609713;Name=ENSE000... 1953399 4 havana exon 165464390 165464586 . + . Parent=transcript:ENST00000511992;Name=ENSE000... In [123]: edf.type.value_counts() Out[123]: exon 1180596 CDS 704604 five_prime_UTR 142387 three_prime_UTR 133938 transcript 96375 gene 42470 processed_transcript 28228 ... Name: type, dtype: int64

isin 메서드는 필터링에 다시 사용됩니다.

그런 다음 빠른 값 카운트는 항목의 대부분이 엑손, 코딩 시퀀스(CDS) 및 번역되지 않은 영역(UTR)임을 보여줍니다.

이것들은 하위 유전자 요소이지만 우리는 주로 유전자 수를 찾고 있습니다. 표시된 대로 42,470개가 있지만 더 알고 싶습니다.

구체적으로 그들의 이름은 무엇이며 어떤 일을 합니까? 이러한 질문에 답하려면 속성 열의 정보를 자세히 살펴보아야 합니다.

 In [127]: ndf = edf[edf.type == 'gene'] In [173]: ndf = ndf.copy() In [133]: ndf.sample(10).attributes.values Out[133]: array(['ID=gene:ENSG00000228611;Name=HNF4GP1;biotype=processed_pseudogene;description=hepatocyte nuclear factor 4 gamma pseudogene 1 [Source:HGNC Symbol%3BAcc:HGNC:35417];gene_id=ENSG00000228611;havana_gene=OTTHUMG00000016986;havana_version=2;logic_name=havana;version=2', 'ID=gene:ENSG00000177189;Name=RPS6KA3;biotype=protein_coding;description=ribosomal protein S6 kinase A3 [Source:HGNC Symbol%3BAcc:HGNC:10432];gene_id=ENSG00000177189;havana_gene=OTTHUMG00000021231;havana_version=5;logic_name=ensembl_havana_gene;version=12', 'ID=gene:ENSG00000231748;Name=RP11-227H15.5;biotype=antisense;gene_id=ENSG00000231748;havana_gene=OTTHUMG00000018373;havana_version=1;logic_name=havana;version=1', 'ID=gene:ENSG00000227426;Name=VN1R33P;biotype=unitary_pseudogene;description=vomeronasal 1 receptor 33 pseudogene [Source:HGNC Symbol%3BAcc:HGNC:37353];gene_id=ENSG00000227426;havana_gene=OTTHUMG00000154474;havana_version=1;logic_name=havana;version=1', 'ID=gene:ENSG00000087250;Name=MT3;biotype=protein_coding;description=metallothionein 3 [Source:HGNC Symbol%3BAcc:HGNC:7408];gene_id=ENSG00000087250;havana_gene=OTTHUMG00000133282;havana_version=3;logic_name=ensembl_havana_gene;version=8', 'ID=gene:ENSG00000177108;Name=ZDHHC22;biotype=protein_coding;description=zinc finger DHHC-type containing 22 [Source:HGNC Symbol%3BAcc:HGNC:20106];gene_id=ENSG00000177108;havana_gene=OTTHUMG00000171575;havana_version=3;logic_name=ensembl_havana_gene;version=5', 'ID=gene:ENSG00000249784;Name=SCARNA22;biotype=scaRNA;description=small Cajal body-specific RNA 22 [Source:HGNC Symbol%3BAcc:HGNC:32580];gene_id=ENSG00000249784;logic_name=ncrna;version=1', 'ID=gene:ENSG00000079101;Name=CLUL1;biotype=protein_coding;description=clusterin like 1 [Source:HGNC Symbol%3BAcc:HGNC:2096];gene_id=ENSG00000079101;havana_gene=OTTHUMG00000178252;havana_version=7;logic_name=ensembl_havana_gene;version=16', 'ID=gene:ENSG00000229224;Name=AC105398.3;biotype=antisense;gene_id=ENSG00000229224;havana_gene=OTTHUMG00000152025;havana_version=1;logic_name=havana;version=1', 'ID=gene:ENSG00000255552;Name=LY6G6E;biotype=protein_coding;description=lymphocyte antigen 6 complex%2C locus G6E (pseudogene) [Source:HGNC Symbol%3BAcc:HGNC:13934];gene_id=ENSG00000255552;havana_gene=OTTHUMG00000166419;havana_version=1;logic_name=ensembl_havana_gene;version=7'], dtype=object)

태그-값 쌍의 세미콜론으로 구분된 목록으로 형식이 지정됩니다. 우리가 가장 관심 있는 정보는 유전자 이름, 유전자 ID, 설명이며 정규식(regex)으로 추출합니다.

 import re RE_GENE_NAME = re.compile(r'Name=(?P<gene_name>.+?);') def extract_gene_name(attributes_str): res = RE_GENE_NAME.search(attributes_str) return res.group('gene_name') ndf['gene_name'] = ndf.attributes.apply(extract_gene_name)

먼저 유전자 이름을 추출합니다.

정규식에서 Name=(?P<gene_name>.+?); , +? 탐욕스럽지 않고 첫 번째 세미콜론에서 검색을 중지하기를 원하기 때문에 + 대신 사용됩니다. 그렇지 않으면 결과는 마지막 세미콜론까지 일치합니다.

또한 정규식은 수천 개의 속성 문자열에 적용할 것이기 때문에 더 나은 성능을 위해 re.search 에서 직접 사용하는 대신 re.compile 로 먼저 컴파일됩니다.

extract_gene_name 은 데이터 프레임 또는 시리즈의 모든 항목에 함수를 적용해야 할 때 사용하는 방법인 pd.apply 에서 사용할 도우미 함수 역할을 합니다.

이 특별한 경우에 ndf.attributes 의 모든 항목에 대한 유전자 이름을 추출하고 gene_name 이라는 새 열의 ndf 에 이름을 다시 추가 gene_name 합니다.

유전자 ID와 설명도 비슷한 방식으로 추출됩니다.

 RE_GENE_ID = re.compile(r'gene_id=(?P<gene_id>ENSG.+?);') def extract_gene_id(attributes_str): res = RE_GENE_ID.search(attributes_str) return res.group('gene_id') ndf['gene_id'] = ndf.attributes.apply(extract_gene_id) RE_DESC = re.compile('description=(?P<desc>.+?);') def extract_description(attributes_str): res = RE_DESC.search(attributes_str) if res is None: return '' else: return res.group('desc') ndf['desc'] = ndf.attributes.apply(extract_description)

RE_GENE_ID 에 대한 정규식은 모든 gene_id 가 ENSG 로 시작해야 한다는 것을 알고 있기 때문에 조금 더 구체적입니다. 여기서 ENS 는 앙상블을 의미하고 G 는 유전자를 의미합니다.

설명이 없는 항목의 경우 빈 문자열을 반환합니다. 모든 것을 추출한 후에는 속성 열을 더 이상 사용하지 않을 .drop 메소드를 사용하여 깔끔하고 깨끗하게 유지하도록 합시다.

 In [224]: ndf.drop('attributes', axis=1, inplace=True) In [225]: ndf.head() Out[225]: seqid source type start end score strand phase gene_id gene_name desc 16 1 havana gene 11869 14409 . + . ENSG00000223972 DDX11L1 DEAD/H-box helicase 11 like 1 [Source:HGNC Sym... 28 1 havana gene 14404 29570 . - . ENSG00000227232 WASH7P WAS protein family homolog 7 pseudogene [Sourc... 71 1 havana gene 52473 53312 . + . ENSG00000268020 OR4G4P olfactory receptor family 4 subfamily G member... 74 1 havana gene 62948 63887 . + . ENSG00000240361 OR4G11P olfactory receptor family 4 subfamily G member... 77 1 ensembl_havana gene 69091 70008 . + . ENSG00000186092 OR4F5 olfactory receptor family 4 subfamily F member...

위의 호출에서 attributes 은 삭제하려는 특정 열을 나타냅니다.

axis=1 은 행 대신 열을 삭제한다는 의미입니다(기본적으로 axis=0 ).

inplace=True 는 지정된 열이 삭제된 새 복사본을 반환하는 대신 DataFrame 자체에서 삭제가 수행됨을 의미합니다.

간단히 .head 를 살펴보면 속성 열이 실제로 사라졌고 세 개의 새로운 열인 gene_name , gene_id 및 desc 가 추가되었음을 알 수 있습니다.

호기심에 모든 gene_id 와 gene_name 이 고유한지 살펴보겠습니다.

 In [232]: ndf.shape Out[232]: (42470, 11) In [233]: ndf.gene_id.unique().shape Out[233]: (42470,) In [234]: ndf.gene_name.unique().shape Out[234]: (42387,)

놀랍게도 유전자 이름의 수가 유전자 ID의 수보다 적어 일부 gene_name이 여러 유전자 ID와 일치해야 함을 나타냅니다. 그들이 무엇인지 알아보자.

 In [243]: count_df = ndf.groupby('gene_name').count().ix[:, 0].sort_values().ix[::-1] In [244]: count_df.head(10) Out[244]: gene_name SCARNA20 7 SCARNA16 6 SCARNA17 5 SCARNA15 4 SCARNA21 4 SCARNA11 4 Clostridiales-1 3 SCARNA4 3 C1QTNF9B-AS1 2 C11orf71 2 Name: seqid, dtype: int64 In [262]: count_df[count_df > 1].shape Out[262]: (63,) In [263]: count_df.shape Out[263]: (42387,) In [264]: count_df[count_df > 1].shape[0] / count_df.shape[0] Out[264]: 0.0014863047632528842

gene_name 값으로 모든 항목을 그룹화한 다음 .count() 를 사용하여 각 그룹의 항목 수를 계산합니다.

ndf.groupby('gene_name').count() 의 출력을 검사하면 모든 열이 각 그룹에 대해 계산되지만 대부분은 동일한 값을 갖습니다.

NA 값은 계산할 때 고려되지 않으므로 첫 번째 열인 seqid 만 계산합니다( NA 값이 없는지 확인하기 위해 .ix[:, 0] 사용).

그런 다음 .sort_values 로 count 값을 정렬하고 .ix[::-1] 로 순서를 반대로 하십시오.

결과적으로 유전자 이름은 최대 7개의 유전자 ID와 공유될 수 있습니다.

 In [255]: ndf[ndf.gene_name == 'SCARNA20'] Out[255]: seqid source type start end score strand phase gene_id gene_name desc 179399 1 ensembl gene 171768070 171768175 . + . ENSG00000253060 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 201037 1 ensembl gene 204727991 204728106 . + . ENSG00000251861 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 349203 11 ensembl gene 8555016 8555146 . + . ENSG00000252778 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 718520 14 ensembl gene 63479272 63479413 . + . ENSG00000252800 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 837233 15 ensembl gene 75121536 75121666 . - . ENSG00000252722 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 1039874 17 ensembl gene 28018770 28018907 . + . ENSG00000251818 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 1108215 17 ensembl gene 60231516 60231646 . - . ENSG00000252577 SCARNA20 small Cajal body-specific RNA 20 [Source:HGNC Symbol%3BAcc:HGNC:32578]

모든 SCARNA20 유전자를 자세히 살펴보면 실제로 모두 다르다는 것을 알 수 있습니다.

그들은 같은 이름을 공유하지만 게놈의 다른 위치에 있습니다.

그러나 그들의 설명은 그들을 구별하는 데 별로 도움이 되지 않는 것 같습니다.

여기서 요점은 유전자 이름이 모든 유전자 ID에 대해 고유하지 않으며 여러 유전자가 공유하는 이름의 약 0.15%를 알고 있다는 것입니다.

전형적인 유전자의 길이는 얼마입니까?

게놈의 불완전성을 이해하려고 할 때 했던 것과 유사하게 ndf 에 length 열을 쉽게 추가할 수 있습니다.

 In [277]: ndf['length'] = ndf.end - ndf.start + 1 In [278]: ndf.length.describe() Out[278]: count 4.247000e+04 mean 3.583348e+04 std 9.683485e+04 min 8.000000e+00 25% 8.840000e+02 50% 5.170500e+03 75% 3.055200e+04 max 2.304997e+06 Name: length, dtype: float64

.describe() 는 길이 값을 기반으로 몇 가지 간단한 통계를 계산합니다.

• 유전자의 평균 길이는 약 36,000 염기입니다.

• 유전자의 중간 길이는 약 5,200 염기 길이입니다.

• 최소 및 최대 유전자 길이는 각각 약 800만 염기 및 230만 염기 길이입니다.

평균이 중앙값보다 훨씬 크기 때문에 길이 분포가 오른쪽으로 치우쳐 있음을 의미합니다. 좀 더 구체적으로 보기 위해 분포를 플로팅해 보겠습니다.

 import matplotlib as plt ndf.length.plot(kind='hist', bins=50, logy=True) plt.show()

Pandas는 Matplotlib에 대한 간단한 인터페이스를 제공하여 DataFrame 또는 시리즈로 플로팅을 매우 편리하게 만듭니다.

이 경우 50개의 빈이 있는 히스토그램 플롯( kind='hist' )을 원하고 y 축이 로그 스케일( logy=True )에 있도록 한다고 말합니다.

히스토그램에서 대부분의 유전자가 첫 번째 빈에 있음을 알 수 있습니다. 그러나 일부 유전자 길이는 2백만 개 이상의 염기가 될 수 있습니다. 그들이 무엇인지 알아보자:

 In [39]: ndf[ndf.length > 2e6].sort_values('length').ix[::-1] Out[39]: seqid source type start end score strand phase gene_name gene_id desc length 2309345 7 ensembl_havana gene 146116002 148420998 . + . CNTNAP2 ENSG00000174469 contactin associated protein-like 2 [Source:HG... 2304997 2422510 9 ensembl_havana gene 8314246 10612723 . - . PTPRD ENSG00000153707 protein tyrosine phosphatase%2C receptor type ... 2298478 2527169 X ensembl_havana gene 31097677 33339441 . - . DMD ENSG00000198947 dystrophin [Source:HGNC Symbol%3BAcc:HGNC:2928] 2241765 440886 11 ensembl_havana gene 83455012 85627922 . - . DLG2 ENSG00000150672 discs large MAGUK scaffold protein 2 [Source:H... 2172911 2323457 8 ensembl_havana gene 2935353 4994972 . - . CSMD1 ENSG00000183117 CUB and Sushi multiple domains 1 [Source:HGNC ... 2059620 1569914 20 ensembl_havana gene 13995369 16053197 . + . MACROD2 ENSG00000172264 MACRO domain containing 2 [Source:HGNC Symbol%... 2057829

보시다시피 가장 긴 유전자는 CNTNAP2로 명명되었으며 이는 contactin related protein-like 2의 약자입니다. 위키피디아 페이지에 따르면,

이 유전자는 7번 염색체의 거의 1.6%를 차지하며 인간 게놈에서 가장 큰 유전자 중 하나입니다.

물론! 우리는 그것을 스스로 확인했습니다. 대조적으로, 가장 작은 유전자는 어떻습니까? 그들은 8개의 베이스만큼 짧을 수 있는 것으로 밝혀졌습니다.

 In [40]: ndf.sort_values('length').head() Out[40]: seqid source type start end score strand phase gene_name gene_id desc length 682278 14 havana gene 22438547 22438554 . + . TRDD1 ENSG00000223997 T cell receptor delta diversity 1 [Source:HGNC... 8 682282 14 havana gene 22439007 22439015 . + . TRDD2 ENSG00000237235 T cell receptor delta diversity 2 [Source:HGNC... 9 2306836 7 havana gene 142786213 142786224 . + . TRBD1 ENSG00000282431 T cell receptor beta diversity 1 [Source:HGNC ... 12 682286 14 havana gene 22449113 22449125 . + . TRDD3 ENSG00000228985 T cell receptor delta diversity 3 [Source:HGNC... 13 1879625 4 havana gene 10238213 10238235 . - . AC006499.9 ENSG00000271544 23

두 극단적인 경우의 길이는 5배 차이(230만 대 8)로, 이는 거대하며 생명체의 다양성 수준을 나타낼 수 있습니다.

단일 유전자는 우리가 탐구하지 않은 대체 스플라이싱(alternative splicing)이라는 과정을 통해 다양한 단백질로 번역될 수 있습니다. 이러한 정보는 GFF3 파일에도 있지만 이 게시물의 범위를 벗어납니다.

염색체 간 유전자 분포

마지막으로 논의하고 싶은 것은 염색체 간의 유전자 분포이며, 이는 두 개의 DataFrame을 결합하는 .merge 방법을 도입하는 예이기도 합니다. 직관적으로 더 긴 염색체는 더 많은 유전자를 수용할 수 있습니다. 그것이 사실인지 봅시다.

 In [53]: ndf = ndf[ndf.seqid.isin(chrs)] In [54]: chr_gene_counts = ndf.groupby('seqid').count().ix[:, 0].sort_values().ix[::-1] Out[54]: chr_gene_counts seqid 1 3902 2 2806 11 2561 19 2412 17 2280 3 2204 6 2154 12 2140 7 2106 5 2002 16 1881 X 1852 4 1751 9 1659 8 1628 10 1600 15 1476 14 1449 22 996 20 965 13 872 18 766 21 541 Y 436 Name: source, dtype: int64

우리는 이전 섹션에서 chrs 변수를 차용하여 어셈블되지 않은 시퀀스를 필터링하는 데 사용했습니다. 출력에 따라 가장 큰 염색체 1이 실제로 가장 많은 유전자를 가지고 있습니다. Y 염색체는 가장 적은 수의 유전자를 가지고 있지만 가장 작은 염색체는 아닙니다.

사실이 아닌 미토콘드리아(MT)에는 유전자가 없는 것 같습니다.

pd.read_csv 에 의해 반환된 첫 번째 DataFrame df 에 대해 조금 더 필터링하면 모든 MT 유전자가 insdc 소스(하바나, ensembl 또는 ensembl_havana의 소스만 고려한 edf 를 생성할 때 필터링됨)에서 가져온 것임을 알 수 있습니다.

 In [134]: df[(df.type == 'gene') & (df.seqid == 'MT')] Out[134]: seqid source type start end score strand phase attributes 2514003 MT insdc gene 648 1601 . + . ID=gene:ENSG00000211459;Name=MT-RNR1;biotype=M... 2514009 MT insdc gene 1671 3229 . + . ID=gene:ENSG00000210082;Name=MT-RNR2;biotype=M... 2514016 MT insdc gene 3307 4262 . + . ID=gene:ENSG00000198888;Name=MT-ND1;biotype=pr... 2514029 MT insdc gene 4470 5511 . + . ID=gene:ENSG00000198763;Name=MT-ND2;biotype=pr... 2514048 MT insdc gene 5904 7445 . + . ID=gene:ENSG00000198804;Name=MT-CO1;biotype=pr... 2514058 MT insdc gene 7586 8269 . + . ID=gene:ENSG00000198712;Name=MT-CO2;biotype=pr... 2514065 MT insdc gene 8366 8572 . + . ID=gene:ENSG00000228253;Name=MT-ATP8;biotype=p... 2514069 MT insdc gene 8527 9207 . + . ID=gene:ENSG00000198899;Name=MT-ATP6;biotype=p... 2514073 MT insdc gene 9207 9990 . + . ID=gene:ENSG00000198938;Name=MT-CO3;biotype=pr... 2514080 MT insdc gene 10059 10404 . + . ID=gene:ENSG00000198840;Name=MT-ND3;biotype=pr... 2514087 MT insdc gene 10470 10766 . + . ID=gene:ENSG00000212907;Name=MT-ND4L;biotype=p... 2514091 MT insdc gene 10760 12137 . + . ID=gene:ENSG00000198886;Name=MT-ND4;biotype=pr... 2514104 MT insdc gene 12337 14148 . + . ID=gene:ENSG00000198786;Name=MT-ND5;biotype=pr... 2514108 MT insdc gene 14149 14673 . - . ID=gene:ENSG00000198695;Name=MT-ND6;biotype=pr... 2514115 MT insdc gene 14747 15887 . + . ID=gene:ENSG00000198727;Name=MT-CYB;biotype=pr...

This example also shows how to combine two conditions during filtering with & ; the logical operator for “or” would be | .

Note that the parentheses around each condition are required, and this part of the syntax in Pandas is different from Python, which would have been literal and and or .

Next, let's borrow the gdf DataFrame from the previous section as a source for the length of each chromosome:

 In [61]: gdf = gdf[gdf.seqid.isin(chrs)] In [62]: gdf.drop(['start', 'end', 'score', 'strand', 'phase' ,'attributes'], axis=1, inplace=True) In [63]: gdf.sort_values('length').ix[::-1] Out[63]: seqid source type length 0 1 GRCh38 chromosome 248956422 1364641 2 GRCh38 chromosome 242193529 1705855 3 GRCh38 chromosome 198295559 1864567 4 GRCh38 chromosome 190214555 1964921 5 GRCh38 chromosome 181538259 2080148 6 GRCh38 chromosome 170805979 2196981 7 GRCh38 chromosome 159345973 2514125 X GRCh38 chromosome 156040895 2321361 8 GRCh38 chromosome 145138636 2416560 9 GRCh38 chromosome 138394717 328938 11 GRCh38 chromosome 135086622 235068 10 GRCh38 chromosome 133797422 483370 12 GRCh38 chromosome 133275309 634486 13 GRCh38 chromosome 114364328 674767 14 GRCh38 chromosome 107043718 767312 15 GRCh38 chromosome 101991189 865053 16 GRCh38 chromosome 90338345 990810 17 GRCh38 chromosome 83257441 1155977 18 GRCh38 chromosome 80373285 1559144 20 GRCh38 chromosome 64444167 1201561 19 GRCh38 chromosome 58617616 2594560 Y GRCh38 chromosome 54106423 1647482 22 GRCh38 chromosome 50818468 1616710 21 GRCh38 chromosome 46709983 2513999 MT GRCh38 chromosome 16569

The columns that are not relevant to the analysis are dropped for clarity.

Yes, .drop can also take a list of columns and drop them altogether in one operation.

Note that the row with a seqid of MT is still there; we will get back to it later. The next operation we will perform is merge the two datasets based on the values of seqid.

 In [73]: cdf = chr_gene_counts.to_frame(name='gene_count').reset_index() In [75]: cdf.head(2) Out[75]: seqid gene_count 0 1 3902 1 2 2806 In [78]: merged = gdf.merge(cdf, on='seqid') In [79]: merged Out[79]: seqid source type length gene_count 0 1 GRCh38 chromosome 248956422 3902 1 10 GRCh38 chromosome 133797422 1600 2 11 GRCh38 chromosome 135086622 2561 3 12 GRCh38 chromosome 133275309 2140 4 13 GRCh38 chromosome 114364328 872 5 14 GRCh38 chromosome 107043718 1449 6 15 GRCh38 chromosome 101991189 1476 7 16 GRCh38 chromosome 90338345 1881 8 17 GRCh38 chromosome 83257441 2280 9 18 GRCh38 chromosome 80373285 766 10 19 GRCh38 chromosome 58617616 2412 11 2 GRCh38 chromosome 242193529 2806 12 20 GRCh38 chromosome 64444167 965 13 21 GRCh38 chromosome 46709983 541 14 22 GRCh38 chromosome 50818468 996 15 3 GRCh38 chromosome 198295559 2204 16 4 GRCh38 chromosome 190214555 1751 17 5 GRCh38 chromosome 181538259 2002 18 6 GRCh38 chromosome 170805979 2154 19 7 GRCh38 chromosome 159345973 2106 20 8 GRCh38 chromosome 145138636 1628 21 9 GRCh38 chromosome 138394717 1659 22 X GRCh38 chromosome 156040895 1852 23 Y GRCh38 chromosome 54106423 436

Since chr_gene_counts is still a Series object, which doesn't support a merge operation, we need to convert it to a DataFrame object first with .to_frame .

.reset_index() converts the original index (ie seqid ) into a new column and resets current index as 0-based incremental numbers.

The output from cdf.head(2) shows what it looks like. Next, we used the .merge method to combine the two DataFrame on the seqid column ( on='seqid' ).

After merging gdf and cdf , the MT entry is still missing. This is because, by default, .merge operates an inner join, while left join, right join, or outer join are available by tuning the how parameter.

Please refer to the documentation for more details.

Later, you may find that there is also a related .join method. .merge and .join are similar yet have different APIs.

According to the official documentation says

The related DataFrame.join method, uses merge internally for the index-on-index and index-on-column(s) joins, but joins on indexes by default rather than trying to join on common columns (the default behavior for merge). If you are joining on index, you may wish to use DataFrame.join to save yourself some typing.

Basically, .merge is more general-purpose and used by .join .

Finally, we are ready to calculate the correlation between chromosome length and gene_count .

 In [81]: merged[['length', 'gene_count']].corr() Out[81]: length gene_count length 1.000000 0.728221 gene_count 0.728221 1.000000

By default .corr calculates the Pearson correlation between all pairs of columns in a dataframe.

But we have only a single pair of columns in this case, and the correlation turns out to be positive – 0.73.

In other words, the larger the chromosome, the more likely it is to have more genes. Let's also plot the two columns after sorting the value pairs by length :

 ax = merged[['length', 'gene_count']].sort_values('length').plot(x='length', y='gene_count',) # add some margin to both ends of x axis xlim = ax.get_xlim() margin = xlim[0] * 0.1 ax.set_xlim([xlim[0] - margin, xlim[1] + margin]) # Label each point on the graph for (s, x, y) in merged[['seqid', 'length', 'gene_count']].sort_values('length').values: ax.text(x, y - 100, str(s)) 

As seen in image above, even though it is a positive correlation overall, it does not hold for all chromosomes. In particular, for Chromosome 17, 16, 15, 14, 13, the correlation is actually negative, meaning the number of genes on the chromosome decreases as the chromosome size increases.

Findings and Future Research

That ends our tutorial on the manipulation of an annotation file for human genome in GFF3 format with the SciPy stack. The tools we've mainly used include IPython, Pandas, and matplotlib. During the tutorial, not only have we learned some of the most common and useful operations in Pandas, we also answered some very interesting questions about our genome. 요약하자면:

1. 인간 게놈의 약 0.37%는 초안이 나온 지 15년이 넘었지만 여전히 불완전합니다.
2. 우리가 사용한 이 특정 GFF3 파일을 기반으로 하는 인간 게놈에는 약 42,000개의 유전자가 있습니다.
3. 유전자의 길이는 수십에서 2백만 개 이상의 염기까지 다양합니다.
4. 유전자는 염색체 사이에 고르게 분포되어 있지 않습니다. 전반적으로 염색체가 클수록 더 많은 유전자를 수용하지만 염색체의 하위 집합에 대해서는 상관 관계가 음수일 수 있습니다.

GFF3 파일은 주석 정보가 매우 풍부하며 표면만 긁었습니다. 추가 탐색에 관심이 있는 경우 다음과 같은 몇 가지 질문을 할 수 있습니다.

1. 유전자는 일반적으로 몇 개의 전사체를 가지고 있습니까? 유전자의 몇 퍼센트가 1개 이상의 전사체를 가지고 있습니까?
2. 성적표에는 일반적으로 몇 개의 isoform이 있습니까?
3. 전사체에는 일반적으로 몇 개의 엑손, CDS 및 UTR이 있습니까? 크기가 어떻게 되나요?
4. 설명 열에 설명된 대로 기능에 따라 유전자를 분류할 수 있습니까?