生信分析人员数据处理脚本实战 | 生信菜鸟团
我前面写到了生信分析人员如何入门linux和perl，后面还会写R和python的总结，但是在这中间我想插入一个脚本实战指南。其实在我前两篇日志里面也重点提到了学习编程语言最重要的就是实战了，也点出了几个关键词。在实际生物信息学数据处理中应用perl和linux，可以借鉴EMBOSS软件套件，fastx-toolkit等基础软件，实现并且模仿该软件的功能。尤其是SMS2/exonerate/里面的一些常见功能,还有DNA2.0 Bioinformatics Toolbox的一些工具。如果你这些名词不懂，请赶快谷歌！！！ 它们做了什么，输入文件是什么，输出文件是什么，你都可以用脚本实现！
你在实现这些功能的时候就必然会融会贯通变量，控制语句，操作符，文件读写等基本编程功能，还会熟悉生物信息学常见数据格式，数据背后的生物学意义。用什么语言都是一样的，千万不要落入语言之争的下乘，也不要纠结于细节。学习是长期过程，尤其是编程这种事情就跟以前的木匠瓦匠一样，是人生技能，跟游戏不一样，不是一时半会就通过了。
如果你英文还不错，推荐看英文的资料，比如下面的DNA2.0 Bioinformatics Toolbox，就可以浏览该网站做了什么，然后自己把同样的文件，对该文件也进行类似的处理。
如果你还是比较熟悉中文，在这里推荐CJ大神总结的一些实际需求，下面都是一些随用随写的脚本，大神都是一句话就搞定了，但是对新手来说，请按部就班的练习！
-1.查看fastq文件读段平均读长、最大读长、最短读长
0.perl命令行粗暴多文件并行处理（每个线程处理一个文件）
1.从fasta文件中提取特定的某个序列(记录)
2.从fasta文件中批量提取序列(记录)
3.Fastq格式转换为fasta格式
4.常规fasta文件去格式为一行id一行seq
5.快速批量提取读段文件的指定序列 (也可用于去格式的fasta文件)
6.读段个数统计
7.fastq质量值格式转换---用于将phred+64数据转为phred+33数据
8.fastq 5'端trimming
9.去除低质量值碱基数量高于N个的reads--用于phred+33数据
10.去除读段序列含未知碱基N超过一定比例的读段
11. 切除读段两端质量值低于给定阈值的部分并丢弃长度低于给定值的记录 新增双端版本 20140831
12.去除低质量值碱基(Q<给定值)所在比例高于(P大于给定值)的读段---用于phred+33数据
13.DNA序列转mRNA序列
14.perl脚本windows和linux间切换
15.window下打印前10行 或者 打印后10行
16.生成批处理用的无后缀file_list
17.fastq中提取特征读段序列
18.fasta格式CDS转为aa（必须有终止密码子）
19.window下面模拟cut命令-提取文本第二列
20.window下合并多个fa文件
21.window下提取匹配到某一模体的fasta序列
22.提取人类基因组注释文件rRNA注释
23.对sort | uniq -c | 的结果频次由高到低排序，有大用
24.fasta格式的DNA序列反向互补
25.一行id一行序列的fa文件格式化为一行id多行序列
26.按fastq文件标签名对读段顺序进行排序---待优化版
27. 替换fq或fa文件记录的id为指定形式
28.提供一个序列名列表逐一替换fasta记录的id29.根据NCBI gene id 即gi号获取GeneBank上的序列
30.根据蛋白gene_id或accession获取其Genebank上的核苷酸序列
31.比较字符串中两个单字符的频次(比如投票0,1或方向F,R)
32.有同学想知道比对上的读段在genome上正反链的分布情况
33.去除全读段所有碱基质量值均低于某个阈值（如20）的读段(支持单端和双端数据)
34.借用pileup文件直接统计测序数据在各染色体上的分布
35.查看sam中uniq mapped比率
36.查看sam中编辑距离分布
37.统计各行平均值或各列平均值
38.将fa文件(尤其基因组文件)分成每个记录一个文件(要求一行id一行seq,见25)
39.批量重命名
40.win下批量去除文件夹内所有文件中的数字
41.统计SAM文件某一标签(BWA结果)
42.提取长度大于1000bp的fa记录
43.批量提取匹配行（正则匹配，强大） ---稍修改即可用于各类模式匹配批量提取，非常强大
44. fasta中有相同id，增加后缀方便blast建库
45. 多个列表文件，比如gene_ids，取样品特异gene_id
46. 直接统计一个序列的GC含量
47. 直接连接几个序列并将小写转换成大写
48. 序列贪吃蛇
49. 随机提取一定比例的fasta 记录或者fastq记录
50. 单行记录随机分组
51. 按照fasta长度排序fasta文件，修改后也可以用于具有某类特征标记的记录排序 (用于大文件，小文件请直接用hash)
52. 双标签区段提取 (使用范围操作符..)
53. 批量从uniprot上下载序列
54. 准备trimmomatic所需的adapter.fa文件
55. 提取fasta文件特定记录的特定区段
56. 获取GO term Level 2的信息
57. 单标签语句块读取 --（方便解析任何行组织文本-fasta fastq blast...）
58. 核酸序列互补配对的子函数
59. 分隔fa文件 fq文件 genebank文件 为数据小文件
60. 序列格式化成每行等长并打印的子函数
61. 从公司返还的注释结果中提取query2gi2GO.table -- for blast3go62. blast2go anno文件转换成blast3go输入文件
63. 提取任意组装结果最长转录本(so-called Unigenes)或者CDS预测结果中最长序列64. 表格类数据，以某一列为keys组成的Group中仅保留其对应某属性（另一列）中值最大的一类 65. 小文件行随机化 66. 打印匹配行及其前'指定数目'行67. 打印匹配行及其后'指定数目'行 68. -n的多个文件区别对待 69. 按照列名提取文件多列 70. 批量提取多个序列多个区段 71. 输出fasta文件每个序列对应的长度 ID\tLength\n72. jar发布前以来外源lib中的jar瘦身73. 依据step长度输出字符串所有后kmer子串74. 基于SAM文件统计ref的每个序列的uniq counts并输出reads的uniq mapped rate统计信息（用于表达谱差异分析 75. 汇总所有counts table并进行无表达补零操作（用于表达谱差异分析76. 保留fastq文件指定长度的读段最优子串77. 输出fasta文件每个记录的A T G C 字数统计78. 合并配对的读段文件fastq 正反读段交错 79. 统计SAM文件 CIGAR的命令 80. fasta文件去除ID行完全重复的记录 81. 合并所有文件的指定列 82. 根据id文件提取第二个文件中多个id匹配行83. 根据某一列的不同值将一个文件分割为多个文件84. 保留高表达或者去除低表达(WGCNA) 85. 表格类数据依据第一列，加和其他所有列，去冗余 86. ghostz比对到nr的表格提取query2gi.table
87. fastqReader
88. Linux下依据 SRA run number下载SRA数据
89. 快速批量统计fq.gz文件行数
90. 格式化mapman结果（mercator）
91. 基因表达量表格做行标准化
92. 基于ID列表提取表格（考虑待提取的表格中有单ID对应多行记录）
93. 文件批量重命名（提供一个重命名列表）
94. perl批量添加fasta文件前缀 （用于多个样本分开组装后合并并用于去冗余等操作）
95. 对表达量表格或者counts表格 依据平均值进行排序
96. 双联表计算卡方值
97. 整理bowtie的比对结果
98. 基于给定列名顺序调整表格列顺序
99. 整理GeneBank文件 （分离地点）
100. 双列文件 整理 为 0-1 交集矩阵
101. 整理bowtie2的比对结果
102. 整理fastqc结果，提取所有样本的读段数
103. 整理STAR比对结果