0

起初，这个线程可能看起来与遗传学有关，但问题实际上是基于 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 文件结构。

谢谢！

4

0 回答 0