生信分析人员数据处理脚本实战



我前面写到了生信分析人员如何入门linux和perl,后面还会写R和python的总结,但是在这中间我想插入一个脚本实战指南。其实在我前两篇日志里面也重点提到了学习编程语言最重要的就是实战了,也点出了几个关键词。在实际生物信息学数据处理中应用perl和linux,可以借鉴EMBOSS软件套件,fastx-toolkit等基础软件,实现并且模仿该软件的功能。尤其是SMS2/exonerate/里面的一些常见功能,还有DNA2.0 Bioinformatics Toolbox的一些工具。如果你这些名词不懂,请赶快谷歌!!! 它们做了什么,输入文件是什么,输出文件是什么,你都可以用脚本实现!

你在实现这些功能的时候就必然会融会贯通变量,控制语句,操作符,文件读写等基本编程功能,还会熟悉生物信息学常见数据格式,数据背后的生物学意义。用什么语言都是一样的,千万不要落入语言之争的下乘,也不要纠结于细节。学习是长期过程,尤其是编程这种事情就跟以前的木匠瓦匠一样,是人生技能,跟游戏不一样,不是一时半会就通过了。

如果你英文还不错,推荐看英文的资料,比如下面的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记录的id

29.根据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 blast3go

62. blast2go anno文件转换成blast3go输入文件




上一篇:残忍的杀戮者皮克顿,49名女性遇害,遗体被喂住
下一篇:没有了