如何解决使用参考 GTF 文件将样本 PLINK tped 文件中的染色体位置转换为基因位置
起初,这个线程可能看起来与遗传有关,但问题实际上是基于 shell 脚本和编程。我是编码新手,所以有人建议我在 SO 中寻求帮助。
我尝试将 NCBI GTF 文件与 PLINK tped 文件相交,目的是将染色体位置切换到 tped 文件(带有标识符、位置和核苷酸的文件)中的基因位置
所以,我做了以下步骤:
1) Got gtf from
https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.39_GRCh38.p13/GCF_000001405.39_GRCh38.p13_genomic.gtf.gz
(由于文件大而复杂,查看文件结构的最佳方式是下载)
2) unzipped and extracted only genes annotation with
awk -F"\t" '$3=="gene"' GCF_000001405.39_GRCh38.p13_genomic.gtf > aaa.gtf
3) converted NC chromosome symbols from NC_xxx to X with
sed -r 's/^NC_[0]*//;s/\.[0-9]\+//;/^N[WT]_|12920/d' aaa.gtf > bbb.gtf
4) converted to bed format (file with positions)
gtf2bed < bbb.gtf > ccc.bed
5) intersected with sample bed file (obtained from Plink tped file)
format of sample bed file:
chr1 183189 183189
chr1 609407 609407
bedtools intersect -wb -a ccc.bed -b sample.bed > ddd.bed
5) removed unwanted symbols (" and ;) from final intersected file with
cat ddd.bed | sed 's/"//g' | sed 's/;//g' > eee.bed
6) printed out only the needed columns from intersected file
cat eee.bed | awk '{print $4,$33,$35,$34,$36,$37}' > fff.bed
最后我有不同列内容的文件:
KLHL17 0 C C
KLHL17 0 A A
PLEKHN1 966179 0 966179 A A
PLEKHN1 966227 0 966227 G G
PERM1 "protein_coding" "C1orf170" gene_synonym 1 976215
PERM1 "protein_coding" "C1orf170" gene_synonym 1 976669
但是想要的结构应该是这样的:
PLEKHN1 966179 0 966179 A A
PLEKHN1 966227 0 966227 G G
PLEKHN1 966272 0 966272 G G
是不是因为 NCBI GTF 文件不一致?它应该成功地与我的样本文件相交,因此我可以将NC位置切换到样本文件中的基因名称和位置并保存tped文件结构。
谢谢!
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。