QCF90이용한 quality control

qcf90 실행 명령어는 다음과 같다.

 

qcf90 --snpfile 50k_for_gs.txt —pedfile pedigree.txt | tee qcf90_01.log

 

명령어 중 “| tee qcf90_01.log”은 실행 중에 화면에 뿌려지는 내용을 qcf90_01.log에 저장하라는 뜻이다. 그러면 실행이 끝나고 천천히 실행 중에 뿌려진 메시지를 자세히 살펴볼 수 있다. 로그 파일을 저장하는 파일로 qcf90.log는 쓰지 말자. qcf90 프로그램이 사용하고 있다.

 

준비한 SNP 유전자형 파일은 다음과 같다.

 

 

준비한 혈통 파일은 다음과 같다.

 

 

다음은 실행 화면이다.

 

 

화면 출력 내용이 저장된 qcf90_01.log의 내용은 다음과 같다.

 

This is qcf90, Version 1.1.0. Yutaka Masuda, Ignacio Aguilar, Andres Legarra, and Ignacy Misztal University of Georgia This is a public binary. The number of genotyped animals is limited to 200000 and the number of markers is limited to 800000. Options from the command-line: --snpfile 50k_for_gs.txt --pedfile pedi2.txt Start saving logs to a file = "qcf90.log" on Feb 5, 2020 14:49:06 Quick checks for pedigree and/or genotypes: pedigree file = pedi2.txt maximum width of ID = 15 * Checking the snp-file format file container = plain file format = text * Checking the first line of marker file number of characters = 51873 number of ID characters = 16 number of leading header = 16 number of markers = 51857 real genotypes (T=yes) = F width of genotype = 1 precision of genotype = 1 * Counting lines of a file file = 50k_for_gs.txt number of lines = 1577 elapsed time = .1 sec. * Preparations hash initialized key character length = 16 internal key length = 4 number of entries = 15318 initial mamory usage = .5 MB Genotypes checks: * Checking the first line of marker file number of characters = 51873 number of ID characters = 16 number of leading header = 16 number of markers = 51857 real genotypes (T=yes) = F width of genotype = 1 precision of genotype = 1 * Checking the format of whole marker file number of lines = 1577 number of errors = 0 elapsed time = 1.3 sec. Genotypes reads: * Genotypes on memory storage format = packed 32-bit integer transposed = yes order = 3242 x 1577 faster read with streaming = no mamory usage = 19.5 MB * Reading whole marker file Reading mode = standard Genotypes / gene content = integer format = (A16,51857I1) elapsed time = 3.8 sec. Pedigree reads and checks: * Loading pedigree to hash table pedigree file = pedi2.txt number of lines = 9146 hash table mamory usage = .5 MB elapsed time = .1 sec. * Checking missing entry in pedigree elapsed time = .0 sec. * Checking sire-dam inconsistency in pedigree elapsed time = .0 sec. * Checking loops in pedigree elapsed time = .0 sec. Genotypes-pedigree consistency: format = general character IDs number of genotypes = 1577 genotyped not in pedigree = 0 (.0 % of genotyped animals) genotyped w/unknwon parent = 849 (53.8 % of genotyped animals) NOTE: More than 50% of genotyped animals have missing pedigree. Please make sure the pedigree file is correct especially when using a XrefID file. Compute statistics for quality-control step1: markers * Allele frequency and call rate for markers: number of markers = 51857 number of genotyped animals = 1577 call rate < 0.99 = 0 (.0 % of markers) call rate < 0.95 = 0 (.0 % of markers) call rate < 0.90 = 0 (.0 % of markers) allele frequency < 0.05 = 9293 (17.9 % of markers) allele frequency < 0.01 = 3602 (6.9 % of markers) monomorphic markers = 1179 (2.3 % of markers) Calculated allele frequency is seved to a file "freqdata.count". elapsed time = 1.0 sec. * Hardy-Weinberg test: This test will use all available marker data. HWE p < 0.05 = 5406 (10.4 % of markers) HWE p < 0.01 = 3510 (6.8 % of markers) HWE p < 0.001 = 2604 (5.0 % of markers) HWE p < 0.0001 = 2283 (4.4 % of markers) HWE p < 0.00001 = 2087 (4.0 % of markers) HWE p < 0.000001 = 1983 (3.8 % of markers) HWE p < .10000E-08 = 1847 (3.6 % of markers) elapsed time = .9 sec. * Summary of unqualified markers: Total number of markers = 51857 QC for call rate: marker call rate < 0.9000 unqualified = 0 (.0 % of markers) QC for minor allele frequency: MAF or (1-MAF) < 0.0500 unqualified = 9293 (17.9 % of markers) QC for monomorphic markers: MAF=1 or MAF=0 unqualified = 0 (.0 % of markers) Some monomorphic markers may be categozied as "low MAF". Total unqualified markers = 9293 (17.9 % of markers) Total qualified markers = 42564 (82.1 % of markers) The unqualified markers will be excluded from the subsequent quality control. elapsed time = .0 sec. Compute statistics for quality-control step2: animals * Call rate for animals: number of effective markers = 42564 number of genotyped animals = 1577 call rate < 0.99 = 0 (.0 % of animals) call rate < 0.95 = 0 (.0 % of animals) call rate < 0.90 = 0 (.0 % of animals) elapsed time = .0 sec. QC for animal call rate: < 0.9000 unqualified = 0 (.0 % of animals) The unqualified animals will be excluded from the subsequent quality control. Log for low call-rate = "Gen_call_rate" * Mendelian conflicts for markers: The sex chromosomes will not be excluded. Number of compared markers = 42564 Comparison mode = Duo QC for marker Mendelian conflicts: > 0.1000 unqualified = 3 (.0 % of all markers) The unqualified markers will be excluded from the subsequent quality control. Number of effective markers = 42561 elapsed time = .1 sec. * Mendelian conflicts for animals: Recalculate conflicts after removing unqualified markers. Animal Mendelian conflicts = "Gen_conflicts" Number of genotyped animals = 1577 Number of effective animals = 1577 Number of compared sires = 393 Number of compared dams = 454 Number of compared trios = 144 Number of effective markers = 42561 Number of conflict pairs = 0 QC for animal Mendelian conflicts: > 0.0100 unqualified = 0 (.0 % of all animals) The animals will be removed from the subsequent quality control. Number of effective animals = 1577 elapsed time = .1 sec. Recalculate allele freqency without unqualified markers and animals * Allele frequency and call rate for markers: Calculated allele frequency is seved to a file "freqdata.count.after.clean". elapsed time = 1.1 sec. Write a status file quality status file = "qcf90.status" Memory usage: hash table = .5 MB integer variables = .5 MB real variables = .4 MB genotypes = 19.5 MB Total memory usage = 20.9 MB Finalizing the program. Done. End logs on Feb 5,2020 14:49:16

 

내용을 하나씩 살펴 보자.

 

Options from the command-line:

--snpfile 50k_for_gs.txt

--pedfile pedi2.txt

 

입력한 SNP 유전자형 파일과 혈통 파일이다.

 

Start saving logs to a file = "qcf90.log" on Feb 5, 2020 14:49:06

 

로그를 qcf90.log에 쓴다고 하나 별 내용이 없다.

 

Quick checks for pedigree and/or genotypes:

pedigree file = pedi2.txt

maximum width of ID = 15

 

혈통 파일의 이름과 혈통 파일에 있는 ID의 최대 길이가 15이다.

 

* Checking the first line of marker file

number of characters = 51873

number of ID characters = 16

number of leading header = 16

number of markers = 51857

real genotypes (T=yes) = F

width of genotype = 1

precision of genotype = 1

 

마커의 수가 51857개이다.

 

* Counting lines of a file

file = 50k_for_gs.txt

number of lines = 1577

elapsed time = .1 sec.

 

SNP 유전자형 파일이름과 1577행이 있다.

 

* Checking sire-dam inconsistency in pedigree

elapsed time = .0 sec.

 

* Checking loops in pedigree

elapsed time = .0 sec.

 

혈통의 오류를 검사하고 있다. 여기서는 나오지 않지만 개체에 나오지 않은 아비와 어미가 있을 경우 에러를 출력한다. 그리고 어떤 개체가 아비나 어미로 동시에 사용되었을 경우 에러를 출력한다. 그리고 자손이 조상으로 나올 경우(loop error) 에러를 출력한다. loop error가 있을 경우 qc_pedigree_loop.doc 파일을 만들고 GraphViz(http://www.graphviz.org)로 읽어 에러를 볼 수 있다.

 

Genotypes-pedigree consistency:

format = general character IDs

number of genotypes = 1577

genotyped not in pedigree = 0 (.0 % of genotyped animals)

genotyped w/unknwon parent = 849 (53.8 % of genotyped animals)

NOTE: More than 50% of genotyped animals have missing pedigree.

Please make sure the pedigree file is correct especially when using a XrefID file.

 

혈통에 없는 개체인데 유전자형을 가지고 있는 개체는 없다. 유전자형을 가진 개체의 50% 이상의 혈통이 missing이다.

 

지금까지는 유전자형 파일과 혈통 파일을 읽은 것이고 qc를 시작한다.

 

Compute statistics for quality-control step1: markers

 

마커에 대한 qc 통계량을 계산한다.

 

* Allele frequency and call rate for markers:

number of markers = 51857

number of genotyped animals = 1577

call rate < 0.99 = 0 (.0 % of markers)

call rate < 0.95 = 0 (.0 % of markers)

call rate < 0.90 = 0 (.0 % of markers)

allele frequency < 0.05 = 9293 (17.9 % of markers)

allele frequency < 0.01 = 3602 (6.9 % of markers)

monomorphic markers = 1179 (2.3 % of markers)

Calculated allele frequency is seved to a file "freqdata.count".

elapsed time = 1.0 sec.

 

1577두의 51857개의 마커에 대한 qc

call rate가 0.99, 0.95 0.90 미만인 마커는 없음. 이것은 imputation을 했기 때문에 당연히 없다. imputation을 하지 않고, finalreport를 변환한 것이라면 어떤 통계가 나온다.

allele frequency가 0.05, 0.01, monomorphic 한 비율이 17.9%, 6.9% 및 2.3%이다. 개량은 차이를 기초로 한다. 이렇게 차이가 없는 마커들은 개량에 의미가 없다는 뜻이다.

그래서 마커별 allele frequency 는 freqdata.count 파일에 있다. freqdata.count는 다음과 같다.

 

 

* Hardy-Weinberg test:

This test will use all available marker data.

HWE p < 0.05 = 5406 (10.4 % of markers)

HWE p < 0.01 = 3510 (6.8 % of markers)

HWE p < 0.001 = 2604 (5.0 % of markers)

HWE p < 0.0001 = 2283 (4.4 % of markers)

HWE p < 0.00001 = 2087 (4.0 % of markers)

HWE p < 0.000001 = 1983 (3.8 % of markers)

HWE p < .10000E-08 = 1847 (3.6 % of markers)

elapsed time = .9 sec.

 

하디 와인버그 평형 검정 결과이다. p < .10000E-08 인 마커의 수가 1847개이다. 이 값이 낮으면 낮을수록 GWAS에서 형질에 많은 영향을 미치는 마커이다.

 

* Summary of unqualified markers:

Total number of markers = 51857

QC for call rate: marker call rate < 0.9000

unqualified = 0 (.0 % of markers)

QC for minor allele frequency: MAF or (1-MAF) < 0.0500

unqualified = 9293 (17.9 % of markers)

QC for monomorphic markers: MAF=1 or MAF=0

unqualified = 0 (.0 % of markers)

Some monomorphic markers may be categozied as "low MAF".

Total unqualified markers = 9293 (17.9 % of markers)

Total qualified markers = 42564 (82.1 % of markers)

The unqualified markers will be excluded from the subsequent quality control.

elapsed time = .0 sec.

 

삭제될 마커의 종류가 개수가 나와 있다. imputation을 해서 삭제될 마커가 별로 없다. MAF 가 5% 미만인 거가 삭제된다. 17.9%의 마커를 삭제하고 개체에 대한 qc를 진행한다.

 

Compute statistics for quality-control step2: animals

 

개체에 대한 qc를 진행한다.

 

* Call rate for animals:

number of effective markers = 42564

number of genotyped animals = 1577

call rate < 0.99 = 0 (.0 % of animals)

call rate < 0.95 = 0 (.0 % of animals)

call rate < 0.90 = 0 (.0 % of animals)

elapsed time = .0 sec.

QC for animal call rate: < 0.9000

unqualified = 0 (.0 % of animals)

The unqualified animals will be excluded from the subsequent quality control.

Log for low call-rate = "Gen_call_rate"

 

역시 imputation을 했기 때문에 call rate가 낮은 개체가 없다. 그러나 지놈스튜디오에서 볼 때 즉 imputation 전에 call rate가 낮은 것(90% 미만 정도, 특정한 기준이 있는 것은 아닌 거 같다.)은 제거하고 연구를 진행하는 것이 정신 건강에 좋다.

개체별 낮은 call rate에 대한 정보는 Gen_call_rate에 나와 있는데, call rate가 낮은 개체가 없으므로 여기서는 빈 파일이다.

 

* Mendelian conflicts for markers:

The sex chromosomes will not be excluded.

Number of compared markers = 42564

Comparison mode = Duo

QC for marker Mendelian conflicts: > 0.1000

unqualified = 3 (.0 % of all markers)

The unqualified markers will be excluded from the subsequent quality control.

Number of effective markers = 42561

elapsed time = .1 sec.

 

마커별 Mendelian conflicts에 대한 정보이다.

 

* Mendelian conflicts for animals:

Recalculate conflicts after removing unqualified markers.

Animal Mendelian conflicts = "Gen_conflicts"

Number of genotyped animals = 1577

Number of effective animals = 1577

Number of compared sires = 393

Number of compared dams = 454

Number of compared trios = 144

Number of effective markers = 42561

Number of conflict pairs = 0

QC for animal Mendelian conflicts: > 0.0100

unqualified = 0 (.0 % of all animals)

The animals will be removed from the subsequent quality control.

Number of effective animals = 1577

elapsed time = .1 sec.

 

개체에 대한 Medelian conflicts이다. 친자감정 결과가 나온다. 이전 pregsf90을 이용한 qc에서 자세히 설명한 바 있다.

 

Recalculate allele freqency without unqualified markers and animals

 

마커, 개체 버릴 건 버리고 다시 계산한다.

 

* Allele frequency and call rate for markers:

Calculated allele frequency is seved to a file "freqdata.count.after.clean".

elapsed time = 1.1 sec.

 

마커에 대한 allele frequency와 call rate 정보를 freqdata.count.after.clean에 저장했다. freqdata.count.after.clean은 다음과 같다.

 

 

Write a status file

quality status file = "qcf90.status"

 

qc 결과를 qcf90.status에 기록했다. qcf90.status는 다음과 같다.

 

 

마커와 개체에 대해서 표시를 했는데 표시에 대한 설명이다.

 

# If you want to create clean files by removing all unqualified markers and animals,

# run qcf90 with the following options.

#

# --statfile qcf90.status --save-clean --snpfile 50k_for_gs.txt --pedfile pedi2.txt

#

# If you want to create clean files by removing markers and animals with secific flags,

# run qcf90 with the above options plus 2 more options. For example, the following options

# will remove the markers with flags 1, 2, and 5 as well as the animals with flag 2.

#

# --statfile qcf90.status --save-clean --snpfile 50k_for_gs.txt --pedfile pedi2.txt --cleanup-marker-flags 1,2,5 --cleanup-animal-flags 2

 

모든 unqualified markers and animals을 제거하고 파일로 저장하기를 바라면 다음과 같은 옵션을 주고 qcf90을 실행시킨다.

 

--statfile qcf90.status --save-clean --snpfile 50k_for_gs.txt —pedfile pedi2.txt

 

특정 플래그를 가진 마커와 개체를 저장하기를 바란다면 다음과 같은 옵션을 주고 qcf90을 실행시킨다.

 

--statfile qcf90.status --save-clean --snpfile 50k_for_gs.txt --pedfile pedi2.txt --cleanup-marker-flags 1,2,5 —cleanup-animal-flags 2

 

마커는 1, 2, 5 플래그를 가진 마커만, 개체는 2 플래그를 가진 개체만 저장한다.

 

마커에 대한 플래그 정보는 다음과 같다.

 

 

개체에 대한 플래그 정보는 다음과 같다.

 

 

쌍둥이인 경우 개체들 사이에 유전자형이 동일할 수 있다. 유전자형이 동일한 지 체크하려면 아래와 같은 옵션을 주고 qcf90을 실행한다.

 

--check-identity

 

qcf90을 통하여 마커와 개체에 대해서 clean하고 다음 단계로 넘어갈 수도 있고, qcf90을 통하여 genotyping이 똑바로 되었는지 검사하고, 즉 혈통과 그에 따른 유전자형이 모두 올바른지 정도만 체크하고 오류를 수정하고 넘어갈 수도 있다. 실제 분석에서는 pregsf90이나 blupf90을 이용하여 처리할 수 있다.