반응형
bedtools 를 이용해서 bed 파일 내 중복되는 영역을 찾아보겠습니다.
bed 파일 외에도 GFF/VCF/BAM 파일 모두 input 으로 가능합니다.
저는 질환과 연관된 병원성 변이 영역이 유전자 영역과 겹치는지 확인해보고 싶었고,
이를 위해 두 가지 .bed 파일을 준비하였습니다.
1. 질환 영역이 포함된 disease.txt
2, 유전자 영역이 포함된 gene.txt
각 .bed 파일의 내용을 간단하게 보시면 아래와 같습니다.
염색체 번호, start, end postion, 질환명 또는 유전자명으로 annotation 되어 있습니다.
[skkwon1048@server test_analysis_blog_221217]$ head disease.txt
chr1 146533376 147883376 1q21.1_microduplication_syndrome
chr1 171808618 173837340 1q24-q25_deletion_syndrome
chr1 209717062 210110985 van_der_Woude_syndrome
chr2 163698001 169702000 2q24.3_deletion_syndrome
chr2 169698001 178002000 2q31.1_duplication_syndrome
chr3 186874563 194994615 3q27.3_deletion_syndrome
chr4 139500001 191154276 4q31-qter_deletion_syndrome
chr4 155598001 164502000 4q32.1-q32.2_triplication_syndrome
chr5 82798001 98202000 5q14.3_deletion_syndrome
chr6 88000001 139000000 6q15-q23_deletions
[skkwon1048@server test_analysis_blog_221217]$ head gene.txt
chr1 55464617 55474465 BSND
chr1 61548295 61922036 NFIA
chr1 55464617 55474465 BSND
chr1 61548295 61922036 NFIA
chr3 189349216 189615062 TP73L
chr3 196466915 196559354 PAK2
chr4 88928820 88998929 PKD2
chr4 111538580 111563117 PITX2
chr5 171288556 171433877 FBXW11
chr5 172659138 172662262 NKX2-5
그럼 이제 bedtools intersect 옵션을 사용해서 겹치는 영역을 보겠습니다.
가장 간단한 default behavior 부터,
bedtools intersect -a disease.txt -b gene.txt
[skkwon1048@server test_analysis_blog_221217]$ bedtools intersect -a disease.txt -b gene.txt
chr3 189349216 189615062 3q27.3_deletion_syndrome
chr6 100836750 100911551 6q15-q23_deletions
chr6 101846954 102516781 6q15-q23_deletions
chr7 74072030 74142672 7q11.23_duplication_syndrome
chr14 57267425 57277184 14q11-q22_deletion_syndrome
chr17 1421283 1466110 17p13.3_Telomeric_duplication_syndrome
chrX 31138746 33367858 Xp21_deletion_syndrome
chrY 25130410 25151610 Partial_chromosome_Y_deletion
chrY 14722879 14823506 Partial_chromosome_Y_deletion
chrY 14542921 15094080 Partial_chromosome_Y_deletion
chrY 25254348 27104406 Partial_chromosome_Y_deletion
[skkwon1048@server test_analysis_blog_221217]$ bedtools intersect -a gene.txt -b disease.txt
chr3 189349216 189615062 TP73L
chr6 100836750 100911551 SIM1
chr6 101846954 102516781 GRIK2
chr7 74072030 74142672 GTF2I
chr14 57267425 57277184 OTX2
chr17 1421283 1466110 PITPNA
chrX 31138746 33367858 DMD
chrY 14542921 15094080 AZFa
chrY 14722879 14823506 Centromere2
chrY 25130410 25151610 BPY2
chrY 25254348 27104406 AZFc-a 옵션이 붙은 파일을 기준으로 overlap 되는 영역을 표시 해줍니다.
[skkwon1048@server test_analysis_blog_221217]$ grep -isrn Partial_chromosome_Y_deletion
disease.txt:39:chrY 200001 28000000 Partial_chromosome_Y_deletion보니까 Partial chromosome Y deletion 질환 영역에 겹치는 부분이 많이 있네요.
질환 영역 자체가 크다 보니까 유전자를 포함하는 부분이 많은 듯 합니다.
[skkwon1048@server test_analysis_blog_221217]$ grep -isrn 25130410
gene.txt:45:chrY 25130410 25151610 BPY2포지션을 grep 으로 찍으서 어떤 유전자가 포함되어 있는지 직접 확인 해볼 수 있을거 같습니다.
그럼 좀더 유용하게 사용될 수 있는 옵션을 실습해보겠습니다.
-wa : reporting the original A feature
-a 로 input 된 정보를 함께 보여줍니다.
-wb : reporting the original B feature
-b 로 input 된 정보를 함께 보여줍니다.
[skkwon1048@server test_analysis_blog_221217]$ bedtools intersect -a disease.txt -b gene.txt -wb
chr3 189349216 189615062 3q27.3_deletion_syndrome chr3 189349216 189615062 TP73L
chr6 100836750 100911551 6q15-q23_deletions chr6 100836750 100911551 SIM1
chr6 101846954 102516781 6q15-q23_deletions chr6 101846954 102516781 GRIK2
chr7 74072030 74142672 7q11.23_duplication_syndrome chr7 74072030 74175022 GTF2I
chr14 57267425 57277184 14q11-q22_deletion_syndrome chr14 57267425 57277184 OTX2
chr17 1421283 1466110 17p13.3_Telomeric_duplication_syndrome chr17 1421283 1466110 PITPNA
chrX 31138746 33367858 Xp21_deletion_syndrome chrX 31138746 33367858 DMD
chrY 25130410 25151610 Partial_chromosome_Y_deletion chrY 25130410 25151610 BPY2
chrY 14722879 14823506 Partial_chromosome_Y_deletion chrY 14722879 14823506 Centromere2
chrY 14542921 15094080 Partial_chromosome_Y_deletion chrY 14542921 15094080 AZFa
chrY 25254348 27104406 Partial_chromosome_Y_deletion chrY 25254348 27104406 AZFc이렇게 (-wa 든) -wb 를 넣어주면 두가지 .bed. 파일 정보를 같이 볼 수 있습니다.
- wo : the amount of overlap between intersecing features
얼마나 overlap 되어 있는지 영역 수치를 알려줍니다.
[skkwon1048@server test_analysis_blog_221217]$ bedtools intersect -a gene.txt -b disease.txt -wo
chr3 189349216 189615062 TP73L chr3 186874563 194994615 3q27.3_deletion_syndrome 265846
chr6 100836750 100911551 SIM1 chr6 88000001 139000000 6q15-q23_deletions 74801
chr6 101846954 102516781 GRIK2 chr6 88000001 139000000 6q15-q23_deletions 669827
chr7 74072030 74175022 GTF2I chr7 72744455 74142672 7q11.23_duplication_syndrome 70642
chr14 57267425 57277184 OTX2 chr14 20000001 58102000 14q11-q22_deletion_syndrome 9759
chr17 1421283 1466110 PITPNA chr17 400001 10702000 17p13.3_Telomeric_duplication_syndrome 44827
chrX 31138746 33367858 DMD chrX 24898001 37602000 Xp21_deletion_syndrome 2229112
chrY 14542921 15094080 AZFa chrY 200001 28000000 Partial_chromosome_Y_deletion 551159
chrY 14722879 14823506 Centromere2 chrY 200001 28000000 Partial_chromosome_Y_deletion 100627
chrY 25130410 25151610 BPY2 chrY 200001 28000000 Partial_chromosome_Y_deletion 21200
chrY 25254348 27104406 AZFc chrY 200001 28000000 Partial_chromosome_Y_deletion 1850058
- c : the number of overlapping features
몇개가 overlapping 되었는지 알려줍니다.
(아까 Partial chromosome Y deletion 는 4개가 나와야겠죠 확인해보시죠)
[skkwon1048@server test_analysis_blog_221217]$ bedtools intersect -a disease.txt -b gene.txt -c
chr1 146533376 147883376 1q21.1_microduplication_syndrome 0
chr1 171808618 173837340 1q24-q25_deletion_syndrome 0
chr1 209717062 210110985 van_der_Woude_syndrome 0
chr2 163698001 169702000 2q24.3_deletion_syndrome 0
chr2 169698001 178002000 2q31.1_duplication_syndrome 0
chr3 186874563 194994615 3q27.3_deletion_syndrome 1
chr4 139500001 191154276 4q31-qter_deletion_syndrome 0
chr4 155598001 164502000 4q32.1-q32.2_triplication_syndrome 0
chr5 82798001 98202000 5q14.3_deletion_syndrome 0
chr6 88000001 139000000 6q15-q23_deletions 2
chr6 130300001 149000000 6q23-q24_deletions 0
chr6 149000001 161000000 6q25_deletion_syndrome 0
chr7 72744455 74142672 7q11.23_duplication_syndrome 1
chr8 77226464 77766239 8q21.11_microdeletion_syndrome 0
chr8 93298001 99002000 8q22.1_duplication_syndrome 0
chr8 93298001 99002000 Nablus_mask-like_facial_syndrome 0
chr9 200001 45000000 9p_deletion_syndrome 0
chr12 860001 34000000 Pallister-Killian_syndrome 0
chr13 40098001 55302000 Axenfeld-Rieger_syndrome_type_2 0
chr14 20000001 58102000 14q11-q22_deletion_syndrome 1
chr14 79300001 106000000 14q31-q32_deletions 0
chr14 89800001 106000000 14q32_deletion_syndrome 0
chr15 33600000 40100000 Idiopathic_generalized_epilepsy 0
chr16 21512062 30199854 16p11.2-p12.2_microdeletion_syndrome 0
chr16 21946524 22467284 16p12.1_microdeletion_syndrome 0
chr16 28833919 30199855 16p11.2_microdeletion_syndrome 0
chr17 400001 10702000 17p13.3_Telomeric_duplication_syndrome 1
chr17 6498001 10702000 17p13.1_deletion_syndrome 0
chr17 16773072 20222149 Potocki-Lupski_syndrome 0
chr17 57598001 61102000 17q23.1-q23.2_duplication_syndrome 0
chr18 2500001 14000000 18p_deletion_syndrome 0
chr18 43500001 61600000 18q_deletion_syndrome 0
chr19 6900001 13900000 19p13.2_microdeletions 0
chr22 37598001 51200000 22q13_duplication_syndrome 0
chrX 24898001 37602000 Xp21_deletion_syndrome 1
chrX 42398001 46402000 Xp11.3_deletion_syndrome 0
chrX 48316644 52759946 Xp11.2_duplication_syndrome 0
chrX 133598001 138002000 Xq26.3_deletion_syndrome 0
chrY 200001 28000000 Partial_chromosome_Y_deletion 4- f : a minimal overlap fraction
최소로 겹치는 영역을 지정해주고 그 수치에 해당하는 영역만 볼 수 있습니다.
저는 10프로로 잡아서 실습해보겠습니다. -f 0.10
[skkwon1048@server test_analysis_blog_221217]$ bedtools intersect -a disease.txt -b gene.txt -f 0.10 -wa -wb
chrX 24898001 37602000 Xp21_deletion_syndrome chrX 31138746 33367858 DMD
-v : reporting the absence of any overlaping features
이 쯤되니 겹치지 않는 부분이 궁금해지는군요.
[skkwon1048@server test_analysis_blog_221217]$ bedtools intersect -a disease.txt -b gene.txt -v
chr1 146533376 147883376 1q21.1_microduplication_syndrome
chr1 171808618 173837340 1q24-q25_deletion_syndrome
chr1 209717062 210110985 van_der_Woude_syndrome
chr2 163698001 169702000 2q24.3_deletion_syndrome
chr2 169698001 178002000 2q31.1_duplication_syndrome
chr4 139500001 191154276 4q31-qter_deletion_syndrome
chr4 155598001 164502000 4q32.1-q32.2_triplication_syndrome
chr5 82798001 98202000 5q14.3_deletion_syndrome
chr6 130300001 149000000 6q23-q24_deletions
chr6 149000001 161000000 6q25_deletion_syndrome
chr8 77226464 77766239 8q21.11_microdeletion_syndrome
chr8 93298001 99002000 8q22.1_duplication_syndrome
chr8 93298001 99002000 Nablus_mask-like_facial_syndrome
chr9 200001 45000000 9p_deletion_syndrome
chr12 860001 34000000 Pallister-Killian_syndrome
chr13 40098001 55302000 Axenfeld-Rieger_syndrome_type_2
chr14 79300001 106000000 14q31-q32_deletions
chr14 89800001 106000000 14q32_deletion_syndrome
chr15 33600000 40100000 Idiopathic_generalized_epilepsy
chr16 21512062 30199854 16p11.2-p12.2_microdeletion_syndrome
chr16 21946524 22467284 16p12.1_microdeletion_syndrome
chr16 28833919 30199855 16p11.2_microdeletion_syndrome
chr17 6498001 10702000 17p13.1_deletion_syndrome
chr17 16773072 20222149 Potocki-Lupski_syndrome
chr17 57598001 61102000 17q23.1-q23.2_duplication_syndrome
chr18 2500001 14000000 18p_deletion_syndrome
chr18 43500001 61600000 18q_deletion_syndrome
chr19 6900001 13900000 19p13.2_microdeletions
chr22 37598001 51200000 22q13_duplication_syndrome
chrX 42398001 46402000 Xp11.3_deletion_syndrome
chrX 48316644 52759946 Xp11.2_duplication_syndrome
chrX 133598001 138002000 Xq26.3_deletion_syndrome겹치지 않는 부분이 왕창 나옵니다. Gene annotation 이 되어 있지 않는 SNP 들을 조사할때 사용될 수 있겠네요.
이 외에도 일대다로 intersect 할 수 도 있습니다.
bedtools intersect -a a.bed -b d1.bed d2.bed d3.bed
아래를 참고하여 주시기 바랍니다.
https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html
intersect — bedtools 2.30.0 documentation
intersect By far, the most common question asked of two sets of genomic features is whether or not any of the features in the two sets “overlap” with one another. This is known as feature intersection. bedtools intersect allows one to screen for overla
bedtools.readthedocs.io
.
반응형
'Bioinformatics' 카테고리의 다른 글
| [Bedtools] BAM to BED 파일 전환 (0) | 2023.01.16 |
|---|---|
| ERBB2 bioinformatics analysis (0) | 2023.01.12 |
| [Bedtools] fasta 파일에서 특정 영역 추출하기 (1) | 2022.10.04 |
| CopywriteR code complete (0) | 2022.09.30 |
| Lecture 8: RNA-sequence Analysis: Expression, Isoforms (0) | 2022.02.22 |
