题目来自生信技能树论坛

image.png

wget ftp://ftp.ensembl.org/pub/release-87/gtf/homo_sapiens/Homo_sapiens.GRCh38.87.chr.gtf.gz
grep -v "#" Homo_sapiens.GRCh38.87.chr.gtf.gz >hg38.gtf.gz
gunzip hg38.gtf.gz

做题之前我们先看看文件的内容是什么

image.png

gtf有9列
1.染色体名

2.注释信息的来源，比如”Genescan”、”Genbank”

等，可以为空，为空用”.”点号代替

3.注释信息的类型，比如Gene、cDNA、mRNA等，或者是SO对应的编号

4、5.开始和结束位置

6.得分，数字，是注释信息可能性的说明，可以是序列相似性比对时的E-values值或者基因预测是的P-values值。”.”表示为空。

7.序列的方向，

+表示正义链, -反义链 , ? 表示未知

8.阅读框：有数字0、1和2。0代表序列的第一个碱基为密码子的第一个碱基，1代表是密码子第二个，2代表第三个。

9.以多个键值对组成的注释信息描述，键与值之间用”=“；

以上图内容简单的理解一下
红线圈出来的范围一个gene信息，范围是11869-14409bp
在这个范围内有两个转录本（黄色部分）：
transcript1 来自基因的11869-14409
这个范围包括3个exon
transcript2来自基因的12010-13670
这个范围包括6个exon

那么本次要做的事情就主要包括一下几个方面吧
1.统计各个染色体基因数量的分布，包括总的基因和各中类型基因数量的分布
2.统计一下基因转录本数量的分布

就这2个吧，依然是分为python，R，shell三种编程方法

1.python
1.统计各个染色体基因数量的分布，包括总的基因和各中类型基因数量的分布
给出我写的脚本

import sys
import collections
args=sys.argv
filename=args[1]
count_gene=collections.OrderedDict()
with open (filename) as fh:
   for line in fh:
      lineL=line.strip().split("\t")
      chr_num="chr"+lineL[0]
      type_info=lineL[2]
      if type_info=="gene":
         All_gene="all"+"_"+type_info
         if chr_num not in count_gene:
            count_gene[chr_num]={}
            count_gene[chr_num][All_gene]=0
         count_gene[chr_num][All_gene]+=1    #计算all gene个数
         descrip_info=lineL[8]
         descrip_info=descrip_info[:-1]      #主要是有的biotype就是最后一部分了，这时候后面得到的b_type包含“;”,所以一开始就把最后的";"去掉
         descrip_infoL=descrip_info.split("; ")
         biotype=descrip_infoL[4]
         biotypeL=biotype.split(" ")
         b_type=biotypeL[1]
         b_type=eval(b_type)
         if b_type not in count_gene[chr_num]:
              count_gene[chr_num][b_type]=0
         count_gene[chr_num][b_type]+=1
for chr_num,countAll in count_gene.items():
   for k,v in countAll.items():
      print(chr_num,k,v)

结果如下

chr1    protein_coding  2052
chr1    TEC 28
chr1    polymorphic_pseudogene  6
chr1    snRNA   219
chr1    antisense   505
chr1    pseudogene  5
chr1    sense_intronic  66
chr1    processed_transcript    31
chr1    all_gene    5194
chr1    sense_overlapping   17
chr1    processed_pseudogene    897
chr1    rRNA    69
chr1    snoRNA  68
chr1    miRNA   126
chr1    lincRNA 576
chr1    IG_V_pseudogene 1
chr1    unprocessed_pseudogene  206
chr1    transcribed_unitary_pseudogene  4
chr1    unitary_pseudogene  26
chr1    scaRNA  14
chr1    misc_RNA    192
chr1    3prime_overlapping_ncRNA    4
chr1    transcribed_processed_pseudogene    33
chr1    bidirectional_promoter_lncRNA   1
chr1    transcribed_unprocessed_pseudogene  48
chr2    protein_coding  1255
chr2    IG_C_gene   1
chr2    TEC 51
chr2    polymorphic_pseudogene  2
chr2    snRNA   157
chr2    antisense   364
chr2    IG_V_gene   46
chr2    pseudogene  1
chr2    transcribed_unitary_pseudogene  3
chr2    processed_transcript    42
chr2    miRNA   105
chr2    sense_overlapping   10
chr2    processed_pseudogene    751
chr2    rRNA    41
chr2    snoRNA  64
chr2    transcribed_unprocessed_pseudogene  39
chr2    all_gene    3971
chr2    lincRNA 577
chr2    IG_J_gene   5
chr2    IG_V_pseudogene 45
chr2    unprocessed_pseudogene  155
chr2    sense_intronic  37
chr2    unitary_pseudogene  7
chr2    scaRNA  7
chr2    misc_RNA    176
chr2    3prime_overlapping_ncRNA    7
chr2    scRNA   1
chr2    transcribed_processed_pseudogene    22
chr3    protein_coding  1076
chr3    TEC 26
chr3    polymorphic_pseudogene  2
chr3    snRNA   133
chr3    antisense   299
chr3    transcribed_unitary_pseudogene  9
chr3    processed_transcript    22
chr3    miRNA   77
chr3    sense_overlapping   7
chr3    processed_pseudogene    616
chr3    rRNA    29
chr3    snoRNA  56
chr3    sRNA    1
chr3    all_gene    3010
chr3    lincRNA 353
chr3    unprocessed_pseudogene  87
chr3    sense_intronic  29
chr3    unitary_pseudogene  9
chr3    scaRNA  2
chr3    misc_RNA    134
chr3    3prime_overlapping_ncRNA    1
chr3    transcribed_processed_pseudogene    19
chr3    transcribed_unprocessed_pseudogene  23
chr4    protein_coding  751
chr4    TEC 37
chr4    polymorphic_pseudogene  1
chr4    snRNA   116
chr4    antisense   207
chr4    pseudogene  1
chr4    sense_intronic  18
chr4    processed_transcript    16
chr4    miRNA   57
chr4    sense_overlapping   9
chr4    processed_pseudogene    561
chr4    ribozyme    1
chr4    rRNA    24
chr4    snoRNA  30
chr4    all_gene    2505
chr4    lincRNA 405
chr4    unprocessed_pseudogene  113
chr4    transcribed_unitary_pseudogene  2
chr4    unitary_pseudogene  2
chr4    scaRNA  1
chr4    misc_RNA    104
chr4    transcribed_processed_pseudogene    31
chr4    bidirectional_promoter_lncRNA   1
chr4    transcribed_unprocessed_pseudogene  17
chr5    protein_coding  876
chr5    all_gene    2868
chr5    lincRNA 490
chr5    rRNA    25
chr5    miRNA   66
chr5    unprocessed_pseudogene  57
chr5    transcribed_unitary_pseudogene  4
chr5    sense_intronic  24
chr5    unitary_pseudogene  4
chr5    snRNA   103
chr5    transcribed_unprocessed_pseudogene  22
chr5    processed_transcript    23
chr5    misc_RNA    119
chr5    TEC 73
chr5    sense_overlapping   10
chr5    transcribed_processed_pseudogene    9
chr5    processed_pseudogene    625
chr5    vaultRNA    1
chr5    antisense   297
chr5    snoRNA  40
chr6    protein_coding  1045
chr6    TEC 36
chr6    polymorphic_pseudogene  3
chr6    snRNA   107
chr6    antisense   238
chr6    non_coding  1
chr6    transcribed_unitary_pseudogene  2
chr6    processed_transcript    8
chr6    miRNA   60
chr6    sense_overlapping   7
chr6    processed_pseudogene    639
chr6    rRNA    27
chr6    snoRNA  39
chr6    all_gene    2863
chr6    lincRNA 360
chr6    unprocessed_pseudogene  105
chr6    sense_intronic  23
chr6    unitary_pseudogene  10
chr6    scaRNA  1
chr6    misc_RNA    105
chr6    3prime_overlapping_ncRNA    2
chr6    transcribed_processed_pseudogene    23
chr6    transcribed_unprocessed_pseudogene  22
chr7    protein_coding  906
chr7    TEC 37
chr7    polymorphic_pseudogene  2
chr7    snRNA   84
chr7    antisense   252
chr7    pseudogene  3
chr7    transcribed_unprocessed_pseudogene  71
chr7    processed_transcript    40
chr7    miRNA   60
chr7    sense_overlapping   12
chr7    processed_pseudogene    577
chr7    TR_D_gene   1
chr7    rRNA    24
chr7    snoRNA  39
chr7    sRNA    1
chr7    all_gene    2867
chr7    lincRNA 286
chr7    unprocessed_pseudogene  189
chr7    sense_intronic  15
chr7    unitary_pseudogene  13
chr7    TR_J_gene   19
chr7    TR_V_pseudogene 17
chr7    TR_C_gene   4
chr7    transcribed_unitary_pseudogene  2
chr7    misc_RNA    143
chr7    TR_V_gene   57
chr7    transcribed_processed_pseudogene    13
chrX    protein_coding  841
chrX    pseudogene  3
chrX    polymorphic_pseudogene  1
chrX    snRNA   84
chrX    misc_RNA    100
chrX    transcribed_unitary_pseudogene  2
chrX    processed_transcript    10
chrX    TEC 14
chrX    sense_overlapping   6
chrX    processed_pseudogene    723
chrX    rRNA    22
chrX    snoRNA  46
chrX    all_gene    2359
chrX    lincRNA 136
chrX    miRNA   105
chrX    unprocessed_pseudogene  108
chrX    sense_intronic  17
chrX    unitary_pseudogene  3
chrX    scaRNA  1
chrX    antisense   101
chrX    3prime_overlapping_ncRNA    1
chrX    transcribed_processed_pseudogene    12
chrX    transcribed_unprocessed_pseudogene  23
chr8    protein_coding  676
chr8    TEC 32
chr8    IG_V_pseudogene 1
chr8    snRNA   84
chr8    antisense   253
chr8    transcribed_unitary_pseudogene  5
chr8    processed_transcript    16
chr8    miRNA   68
chr8    sense_overlapping   8
chr8    processed_pseudogene    468
chr8    rRNA    29
chr8    polymorphic_pseudogene  1
chr8    all_gene    2353
chr8    lincRNA 412
chr8    unprocessed_pseudogene  114
chr8    sense_intronic  45
chr8    unitary_pseudogene  4
chr8    misc_RNA    82
chr8    3prime_overlapping_ncRNA    1
chr8    transcribed_processed_pseudogene    14
chr8    snoRNA  33
chr8    transcribed_unprocessed_pseudogene  7
chr9    protein_coding  781
chr9    TEC 18
chr9    polymorphic_pseudogene  2
chr9    snRNA   65
chr9    misc_RNA    96
chr9    IG_C_pseudogene 1
chr9    pseudogene  1
chr9    transcribed_unprocessed_pseudogene  28
chr9    processed_transcript    17
chr9    miRNA   74
chr9    sense_overlapping   9
chr9    processed_pseudogene    461
chr9    ribozyme    2
chr9    rRNA    19
chr9    snoRNA  31
chr9    all_gene    2242
chr9    lincRNA 275
chr9    IG_V_pseudogene 4
chr9    unprocessed_pseudogene  132
chr9    sense_intronic  22
chr9    unitary_pseudogene  5
chr9    TR_V_pseudogene 4
chr9    scaRNA  1
chr9    transcribed_unitary_pseudogene  2
chr9    antisense   165
chr9    3prime_overlapping_ncRNA    1
chr9    transcribed_processed_pseudogene    23
chr9    TR_V_gene   3
chr11   protein_coding  1276
chr11   TEC 76
chr11   polymorphic_pseudogene  22
chr11   snRNA   71
chr11   antisense   326
chr11   transcribed_unprocessed_pseudogene  30
chr11   processed_transcript    20
chr11   pseudogene  2
chr11   sense_overlapping   17
chr11   all_gene    3235
chr11   processed_pseudogene    473
chr11   ribozyme    1
chr11   rRNA    24
chr11   snoRNA  52
chr11   sRNA    2
chr11   miRNA   81
chr11   lincRNA 305
chr11   unprocessed_pseudogene  279
chr11   sense_intronic  40
chr11   unitary_pseudogene  20
chr11   scaRNA  3
chr11   transcribed_unitary_pseudogene  3
chr11   misc_RNA    97
chr11   3prime_overlapping_ncRNA    1
chr11   transcribed_processed_pseudogene    14
chr10   protein_coding  732
chr10   pseudogene  1
chr10   IG_V_pseudogene 1
chr10   snRNA   86
chr10   misc_RNA    89
chr10   transcribed_unprocessed_pseudogene  39
chr10   processed_transcript    29
chr10   TEC 31
chr10   sense_overlapping   9
chr10   processed_pseudogene    422
chr10   ribozyme    2
chr10   rRNA    29
chr10   snoRNA  28
chr10   polymorphic_pseudogene  1
chr10   all_gene    2204
chr10   lincRNA 280
chr10   miRNA   61
chr10   unprocessed_pseudogene  85
chr10   sense_intronic  33
chr10   unitary_pseudogene  5
chr10   transcribed_unitary_pseudogene  2
chr10   antisense   226
chr10   transcribed_processed_pseudogene    13
chr12   protein_coding  1034
chr12   TEC 99
chr12   snRNA   99
chr12   antisense   304
chr12   non_coding  1
chr12   transcribed_unprocessed_pseudogene  18
chr12   processed_transcript    22
chr12   miRNA   62
chr12   sense_overlapping   5
chr12   processed_pseudogene    490
chr12   rRNA    27
chr12   snoRNA  37
chr12   all_gene    2940
chr12   lincRNA 438
chr12   unprocessed_pseudogene  69
chr12   sense_intronic  75
chr12   unitary_pseudogene  6
chr12   scaRNA  2
chr12   transcribed_unitary_pseudogene  5
chr12   misc_RNA    115
chr12   3prime_overlapping_ncRNA    2
chr12   macro_lncRNA    1
chr12   transcribed_processed_pseudogene    29
chr13   protein_coding  327
chr13   unprocessed_pseudogene  52
chr13   all_gene    1304
chr13   lincRNA 241
chr13   snRNA   41
chr13   miRNA   31
chr13   antisense   105
chr13   sense_intronic  39
chr13   unitary_pseudogene  1
chr13   transcribed_unitary_pseudogene  1
chr13   transcribed_unprocessed_pseudogene  15
chr13   processed_transcript    12
chr13   misc_RNA    75
chr13   TEC 29
chr13   snoRNA  17
chr13   processed_pseudogene    295
chr13   rRNA    15
chr13   transcribed_processed_pseudogene    8
chr14   protein_coding  623
chr14   IG_C_gene   9
chr14   TEC 20
chr14   polymorphic_pseudogene  3
chr14   TR_J_pseudogene 4
chr14   snRNA   64
chr14   IG_pseudogene   1
chr14   antisense   179
chr14   IG_V_gene   49
chr14   IG_C_pseudogene 2
chr14   non_coding  1
chr14   pseudogene  1
chr14   transcribed_unprocessed_pseudogene  14
chr14   processed_transcript    21
chr14   miRNA   91
chr14   ribozyme    2
chr14   IG_J_pseudogene 3
chr14   processed_pseudogene    324
chr14   TR_D_gene   3
chr14   rRNA    10
chr14   snoRNA  72
chr14   bidirectional_promoter_lncRNA   1
chr14   all_gene    2224
chr14   lincRNA 297
chr14   IG_J_gene   6
chr14   IG_V_pseudogene 84
chr14   unprocessed_pseudogene  47
chr14   sense_intronic  21
chr14   unitary_pseudogene  2
chr14   TR_J_gene   60
chr14   TR_V_pseudogene 9
chr14   TR_C_gene   2
chr14   scaRNA  2
chr14   transcribed_unitary_pseudogene  2
chr14   misc_RNA    79
chr14   IG_D_gene   27
chr14   transcribed_processed_pseudogene    30
chr14   sense_overlapping   11
chr14   TR_V_gene   48
chr15   protein_coding  597
chr15   TEC 47
chr15   IG_V_pseudogene 6
chr15   snRNA   64
chr15   misc_RNA    93
chr15   IG_V_gene   6
chr15   pseudogene  1
chr15   sense_intronic  79
chr15   processed_transcript    25
chr15   all_gene    2152
chr15   sense_overlapping   8
chr15   processed_pseudogene    287
chr15   rRNA    12
chr15   snoRNA  114
chr15   miRNA   57
chr15   lincRNA 300
chr15   unprocessed_pseudogene  118
chr15   IG_D_gene   10
chr15   unitary_pseudogene  1
chr15   scaRNA  3
chr15   transcribed_unitary_pseudogene  3
chr15   antisense   225
chr15   3prime_overlapping_ncRNA    2
chr15   transcribed_processed_pseudogene    29
chr15   transcribed_unprocessed_pseudogene  65
chr16   protein_coding  866
chr16   TEC 136
chr16   IG_V_pseudogene 8
chr16   snRNA   52
chr16   antisense   307
chr16   IG_V_gene   6
chr16   pseudogene  1
chr16   transcribed_unprocessed_pseudogene  53
chr16   processed_transcript    31
chr16   all_gene    2511
chr16   sense_overlapping   12
chr16   processed_pseudogene    252
chr16   rRNA    32
chr16   snoRNA  33
chr16   polymorphic_pseudogene  1
chr16   miRNA   69
chr16   lincRNA 360
chr16   unprocessed_pseudogene  112
chr16   sense_intronic  83
chr16   unitary_pseudogene  12
chr16   scaRNA  1
chr16   transcribed_unitary_pseudogene  2
chr16   misc_RNA    51
chr16   3prime_overlapping_ncRNA    5
chr16   bidirectional_promoter_lncRNA   1
chr16   transcribed_processed_pseudogene    25
chr17   transcribed_unprocessed_pseudogene  64
chr17   protein_coding  1197
chr17   unprocessed_pseudogene  105
chr17   TEC 99
chr17   lincRNA 338
chr17   snRNA   86
chr17   miRNA   75
chr17   misc_RNA    99
chr17   transcribed_unitary_pseudogene  3
chr17   unitary_pseudogene  6
chr17   scaRNA  5
chr17   sense_intronic  69
chr17   processed_transcript    48
chr17   antisense   376
chr17   all_gene    2995
chr17   sense_overlapping   3
chr17   transcribed_processed_pseudogene    43
chr17   processed_pseudogene    310
chr17   rRNA    17
chr17   snoRNA  52
chr18   protein_coding  270
chr18   TEC 50
chr18   snRNA   46
chr18   misc_RNA    41
chr18   IG_C_pseudogene 1
chr18   transcribed_unitary_pseudogene  1
chr18   processed_transcript    11
chr18   miRNA   30
chr18   sense_overlapping   1
chr18   processed_pseudogene    207
chr18   rRNA    13
chr18   snoRNA  17
chr18   sRNA    1
chr18   all_gene    1170
chr18   lincRNA 246
chr18   unprocessed_pseudogene  14
chr18   sense_intronic  49
chr18   unitary_pseudogene  3
chr18   scaRNA  2
chr18   antisense   146
chr18   transcribed_processed_pseudogene    14
chr18   transcribed_unprocessed_pseudogene  7
chr20   protein_coding  544
chr20   all_gene    1386
chr20   lincRNA 214
chr20   snRNA   46
chr20   miRNA   39
chr20   unprocessed_pseudogene  32
chr20   sense_intronic  20
chr20   unitary_pseudogene  5
chr20   transcribed_unitary_pseudogene  3
chr20   scaRNA  1
chr20   transcribed_unprocessed_pseudogene  8
chr20   processed_transcript    9
chr20   misc_RNA    68
chr20   TEC 10
chr20   antisense   138
chr20   transcribed_processed_pseudogene    3
chr20   processed_pseudogene    198
chr20   sense_overlapping   3
chr20   rRNA    20
chr20   snoRNA  25
chr19   protein_coding  1470
chr19   TEC 74
chr19   lincRNA 236
chr19   snRNA   29
chr19   miRNA   115
chr19   unprocessed_pseudogene  159
chr19   polymorphic_pseudogene  2
chr19   sense_intronic  55
chr19   misc_RNA    61
chr19   unitary_pseudogene  7
chr19   transcribed_unprocessed_pseudogene  53
chr19   processed_transcript    36
chr19   antisense   293
chr19   3prime_overlapping_ncRNA    2
chr19   all_gene    2926
chr19   sense_overlapping   6
chr19   snoRNA  22
chr19   processed_pseudogene    263
chr19   rRNA    13
chr19   transcribed_processed_pseudogene    30
chrY    protein_coding  53
chrY    unprocessed_pseudogene  216
chrY    all_gene    523
chrY    lincRNA 52
chrY    snRNA   17
chrY    antisense   5
chrY    sense_intronic  1
chrY    transcribed_unprocessed_pseudogene  27
chrY    processed_transcript    1
chrY    misc_RNA    7
chrY    snoRNA  3
chrY    processed_pseudogene    130
chrY    rRNA    7
chrY    transcribed_processed_pseudogene    4
chr22   protein_coding  438
chr22   IG_C_gene   4
chr22   TEC 9
chr22   IG_V_pseudogene 43
chr22   snRNA   26
chr22   antisense   134
chr22   IG_V_gene   37
chr22   IG_C_pseudogene 5
chr22   transcribed_unprocessed_pseudogene  36
chr22   processed_transcript    18
chr22   miRNA   32
chr22   sense_overlapping   9
chr22   processed_pseudogene    156
chr22   rRNA    6
chr22   snoRNA  11
chr22   polymorphic_pseudogene  2
chr22   all_gene    1318
chr22   lincRNA 165
chr22   IG_J_gene   7
chr22   unprocessed_pseudogene  75
chr22   sense_intronic  31
chr22   unitary_pseudogene  3
chr22   scaRNA  3
chr22   misc_RNA    62
chr22   transcribed_processed_pseudogene    6
chr21   protein_coding  233
chr21   all_gene    835
chr21   lincRNA 191
chr21   snRNA   21
chr21   miRNA   26
chr21   unprocessed_pseudogene  32
chr21   IG_V_gene   1
chr21   sense_intronic  17
chr21   pseudogene  1
chr21   transcribed_unprocessed_pseudogene  6
chr21   processed_transcript    7
chr21   misc_RNA    24
chr21   3prime_overlapping_ncRNA    1
chr21   TEC 16
chr21   antisense   81
chr21   transcribed_processed_pseudogene    3
chr21   processed_pseudogene    143
chr21   sense_overlapping   8
chr21   rRNA    9
chr21   snoRNA  15
chrMT   Mt_rRNA 2
chrMT   protein_coding  13
chrMT   Mt_tRNA 22
chrMT   all_gene    37

可视化一下吧
因为all_gene值太高，所以把它单独拿出来

grep "all_gene" count_all.txt >count_gene.txt

library(ggplot2)
data <- read.table("C:/Users/pc/Desktop/count_gene.txt")
chr_name <- data$V1
type <- data$V2 
number <- data$V3
p <- ggplot(data,aes(x=chr_name,y=number))+theme_bw()+theme(axis.text.x=element_text(angle=90,size =8,vjust = 0.5,hjust = 0.5),panel.grid = element_blank())+geom_bar(aes(fill=chr_name),stat = "identity",position="dodge",width=0.8)

image.png

所有类型基因的可视化如下
看一下一号染色体：

grep "chr1  " count_else.txt >count_else.chr1.txt

image.png

编码蛋白的基因数量还是最多的

最后show一下怎么把多个染色体的图堆在一起，用到ggplot2中的facet

 grep "chr2" count_else.txt >count_else.chr2_all.txt
#取出2，20，22号染色体

image.png

这里只是想展示一下这种画图方法，但是如果把所有染色体都放上就压缩了，图看不清楚。

2.统计一下基因转录本数量的分布
脚本如下：

import sys
import collections
args=sys.argv
filename=args[1]
count_transcript=collections.OrderedDict()
with open (filename) as fh:
   for line in fh:
      lineL=line.strip().split("\t")
      chr_num="chr"+lineL[0]
      type_info=lineL[2]
      if type_info=="transcript":
         if chr_num not in count_transcript:
            count_transcript[chr_num]={}
         descrip_info=lineL[8]
         descrip_infoL=descrip_info.split("; ")
         gene_ID=descrip_infoL[0]
         gene_IDL=gene_ID.split(" ")
         gene_ID=gene_IDL[1]
         gene_ID=eval(gene_ID)
         if gene_ID not in count_transcript[chr_num]:
              count_transcript[chr_num][gene_ID]=0
         count_transcript[chr_num][gene_ID]+=1
for chr_num,countAll in count_transcript.items():
   for k,v in countAll.items():
      print(chr_num,k,v)

2.R
（略）
3.shell
（略）

life is short ,还是把python练习好，R和shell打辅助吧