task4c-研究最热门的基因是什么

作业链接

作业题目

gene2pubmed.gz 约50M文件里面的信息太丰富了，有1333万行信息，仅仅是人类就有159万行的文献，涉及到3万9千的基因数量，绝大部分基因都是如过眼云烟，很少人去研究它。
我们的TP53能拔得头彩也是不容易，但它也有自己的发展规律，希望大家可以更细致去探索 ftp://ftp.ncbi.nlm.nih.gov//gene 里面的文件。
比如这样的top 100的基因词云，其实可以做出来最近30年的变化规律，只需要你去找到文献的时间年份信息，进行拆分，每个年份独立统计绘图即可。

entrez ID 和年份信息提取

基因entrez ID和pubmed id的对应信息 下载:ftp://ftp.ncbi.nlm.nih.gov//gene/DATA/gene2pubmed.gz

$ zcat gene2pubmed.gz |awk 'NR==1 || $1=="9606"' |wc -l
1653326 #目前pubmed中人类的文献有1653325篇文献

pubmed id和文章发表时间的对应信息 下载链接:https://ftp.ncbi.nlm.nih.gov/pubmed/baseline/
把该ftp下的xml.gz 和xml.gz.md5两种文件全部下载下来了，一共大概30G，2124个文件
*xml.gz文件1062个
​

把md5文件合并成一个md5.txt文件

cat *md5 |perl -nle '{/MD5\((.*?)\)=\s+(\w+$)/;print "$2\t$1"}' >md5.txt

$ head md5.txt 
02960ff0fa7de8ca936d93e66c4d1474	pubmed21n0001.xml.gz
5f6439782c8344314f53aacfa9fbbcee	pubmed21n0002.xml.gz
14cd080a2913d73460e0d9be36a40c38	pubmed21n0003.xml.gz
32a203915334ab883b74212fb776f687	pubmed21n0004.xml.gz
6e410e4d62ac622d2a4a24073cbe15c3	pubmed21n0005.xml.gz
185427245e984126fed1bd1239aa782c	pubmed21n0006.xml.gz
144eb0f8df231faec6f27e8f1d07bc70	pubmed21n0007.xml.gz
3ad439edfcd3987e1298baf84a41f6c0	pubmed21n0008.xml.gz
97267082bf3fed1d48e78559a6d12b65	pubmed21n0009.xml.gz
b254ae36145761118eb10a17c9d13273	pubmed21n0010.xml.gz

md5值检验-确保下载下来的数据的完整性

$ md5sum -c md5.txt 
pubmed21n0001.xml.gz: 确定
pubmed21n0002.xml.gz: 确定
pubmed21n0003.xml.gz: 确定
pubmed21n0004.xml.gz: 确定
pubmed21n0005.xml.gz: 确定
pubmed21n0006.xml.gz: 确定
pubmed21n0007.xml.gz: 确定
pubmed21n0008.xml.gz: 确定
pubmed21n0009.xml.gz: 确定
......

xml数据格式
一大堆xml.gz的数据，接下来的问题就是如何 快速便捷的提取pubmed id 和 文章发表年份的信息。

$ ls data/*gz |head
data/pubmed21n0001.xml.gz
data/pubmed21n0002.xml.gz
data/pubmed21n0003.xml.gz
data/pubmed21n0004.xml.gz
data/pubmed21n0005.xml.gz
data/pubmed21n0006.xml.gz
data/pubmed21n0007.xml.gz
data/pubmed21n0008.xml.gz
data/pubmed21n0009.xml.gz
data/pubmed21n0010.xml.gz

解析脚本-xmlParse.pl

$/="<PubmedArticle>";

open F,"gzip -dc $ARGV[0]|";
while(<F>){
	next if $.==1;
	chomp;
	$pmid="";$year="";

	if(/\<PMID.*?\>(\d+)\<\/PMID\>/){$pmid=$1;}

	if(/<PubDate>.*?<Year>(\d+)<\/Year>.*?<\/PubDate>/s){
		$year=$1;
	}elsif(/<DateCompleted>.*?<Year>(\d+)<\/Year>.*?<\/DateCompleted>/s){
		$year=$1;
	}
		
	print "$pmid\t$year\n";
}
close F;

#批量解析并存储
ls data/*gz |while read id;do(perl xmlParse.pl $id >>xml_20211001.txt);done

提取的标签如下
pubmed id

<ArticleId IdType="pubmed">1</ArticleId>

time
有PubDate这个标签就用这个；没有就用DateCompleted这个标签；再没有就空着

解析好的数据:
与gene2pubmed.gz做了一个merge；
只保留了有时间年份的人类 gene_id 和 年份的信息
时间跨度从1902-2021；涉及的entrenz id总数： 40958
链接：https://pan.baidu.com/s/1ovDiraQvWjd6axin8GmDyA 提取码：owkx 大小:16M

$ tail Homo_GeneID_year.txt 
122152361	2016
122152362	2016
122198330	1998
122394733	2003
122394733	2004
122394733	2005
122394733	2009

$ head Homo_GeneID_year.txt 
GeneID	year
1	1989
1	1986
1	1987
1	1996
1	2002
1	2004
1	2004
1	2004
1	2004

Homo_GeneID_PubMed_ID_year.txt #人类 gene_id,pubmed id; 年份信息
链接：https://pan.baidu.com/s/1pixrnK1WsNa6RFZUF_PTUw
提取码：pg6m
大小:38M
比上个文件多两列信息

$ head -5 Homo_GeneID_PubMed_ID_year.txt 
#tax_id	GeneID	PubMed_ID	year
9606	1	2591067	1989
9606	1	3458201	1986
9606	1	3610142	1987
9606	1	8889549	1996

下一步就是gene_id的转换以及可视化了
数据过滤

library(dplyr)
gene_id_year = read.table('Homo_GeneID_year.txt', header = T)

gene_id_last30year = dplyr::filter(gene_id_year, year>=1990) %>% table %>% as.data.frame
stat_gene = dplyr::filter(gene_id_year, year>=1990) %>% table %>% as.data.frame %>% dplyr::filter(.,Freq>0)

gene_id转换

library(org.Hs.eg.db)
ids = toTable(org.Hs.egSYMBOL)
colnames(ids) = c('GeneID', "symbol")
data_merge = dplyr::left_join(stat_gene, ids ,by='GeneID')
data_merge = na.omit(data_merge)

TP53各年份的文献数量展示

library(dplyr)

tp53_data = dplyr::filter(data_merge,symbol=='TP53',year !='2021')%>%dplyr::select(.,symbol,year,Freq)

ggplot(tp53_data, aes(x=year, y=Freq )) +
  geom_bar(stat="identity", fill=rainbow(length(tp53_data$year)),alpha=.6, width=.4) +
  coord_flip() +
  theme_bw()+labs(title = "TP53")