GEO数据挖掘

复现此篇文章: https://mp.weixin.qq.com/s/5Bdx7YyAk8InYs2MFIUUVQ

下载相应的R包

options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/") # 设置镜像
BiocManager::install('stringi')
BiocManager::install('GEOquery')
BiocManager::install('limma')
BiocManager::install('ggfortify')
BiocManager::install('ggplot2')
BiocManager::install('pheatmap')
BiocManager::install('ggstatsplot')
BiocManager::install('VennDiagram')
BiocManager::install('clusterProfiler')
BiocManager::install('org.Hs.eg.db')

浏览器下载gtf格式文件：没录全屏，最后点击链接后，浏览器就开始下载文件了。

gtf.gif

---

gtf格式：GTF基因注释文件详解

gunzip .gtf.zip文件
cat .gtf文件
tail  .gtf文件  

通过以上步骤复制粘贴出前后的两条注释，做初步了解：

1   havana  gene    11869   14409   .   +   .   gene_id "ENSG00000223972"; gene_version "5"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene";

KI270713.1  ensembl CDS 35407   35913   .   +   0   gene_id "ENSG00000268674"; gene_version "2"; transcript_id "ENST00000601199"; transcript_version "2"; exon_number "1"; gene_name "AC213203.1"; gene_source "ensembl"; gene_biotype "protein_coding"; transcript_name "AC213203.1-201"; transcript_source "ensembl"; transcript_biotype "protein_coding"; protein_id "ENSP00000473163"; protein_version "

1. seq_id：序列的编号，一般为chr或者scanfold编号；
2. source: 注释的来源，一般为数据库或者注释的机构，如果未知，则用点“.”代替；
3. type: 注释信息的类型，比如Gene、cDNA、mRNA、CDS等
4. start:该基因或转录本在参考序列上的起始位置；
5. end: 该基因或转录本在参考序列上的终止位置；
6. score: 得分，数字，是注释信息可能性的说明，可以是序列相似性比对时的E-values值或者基因预测是的P-values值，“.”表示为空；
7. strand: 该基因或转录本位于参考序列的正链(+)或负链(-)上;
8. phase: 仅对注释类型为“CDS”有效，表示起始编码的位置，有效值为0、1、2(对于编码蛋白质的CDS来说，本列指定下一个密码子开始的位置。每3个核苷酸翻译一个氨基酸，从0开始，CDS的起始位置，除以3，余数就是这个值，，表示到达下一个密码子需要跳过的碱基个数。该编码区第一个密码子的位置，取值0,1,2。0表示该编码框的第一个密码子第一个碱基位于其5'末端；1表示该编码框的第一个密码子的第一个碱基位于该编码区外；2表示该编码框的第一个密码子的第一、二个碱基位于该编码区外；如果Feature为CDS时，必须指明具体值。)；
9. attributes:一个包含众多属性的列表，格式为“标签＝值”（tag=value），标签与值之间以空格分开，且每个特征之后都要有分号；（包括最后一个特征），其内容必须包括gene_id和transcript_id。以多个键值对组成的注释信息描述，键与值之间用“=”，不同的键值用“；

要得到基因与基因类型的对应关系即需要取出第三列和 gene_name这一列。

cat Homo_sapiens.GRCh38.96.gtf >hg38
awk '{print $3}' hg38 >type1  #输出第三列内容到type1里
sed '/^\s*$/d' type1 >type2  删除空行并输出内容到type2里
 uniq type2 >type3   #去掉重复行
awk '{print NR $0}' type3 | tail -1 #统计行数

awk命令：awk '{pattern + action}' 或者 awk 'pattern {action}'

---

• 作用：把文件逐行的读入，以空格为默认分隔符将每行切片，切开的部分再进行各种分析处理。

• pattern 表示 AWK 在数据中查找的内容，而 action 是在找到匹配内容时所执行的一系列命令。花括号 ({}) 不需要在程序中始终出现，但它们用于根据特定的模式对一系列指令进行分组。

• 变量NF : field 在文本文件表示一个竖行
  NR：record 在文本文件表示一个横行

awk '{print $0}'  tmp.txt         # 打印出tmp.txt的全部内容  $0:全部内容 $1:第一列 $n：第n列
awk'$3==1986{print $0}'  tmp.txt  #第三列是1986则打印出所有内容
awk '{if($0!="") print}' tmp.txt  #删除空行
awk '{if(!NF || /^#/){next}}1' file1 实现对文件里面的空行和#开头的行进行跳过操作，并输出结果。

sed命令：

附：

• 推荐b站awk系列视频
• sed文章

awk文章