perl 一行命令 2021-09-02

yangkfor

已于 2022-07-21 17:25:34 修改

阅读量7.7k

点赞数 2

分类专栏：生物信息学文章标签： perl

于 2021-09-02 13:35:56 首次发布

原文链接：https://user.qzone.qq.com/120509419/blog/1405613131?_t_=0.33067108750695295

版权

生物信息学专栏收录该内容

8 篇文章 3 订阅

订阅专栏

生物信息入门常用

perl ‘one-liner’–by CJ chen
convert -delay 10 -resample 10% ls|perl -lane '/CJ\d{4}\d{4}(\d{2})(\d{2})/;print if $1>8 and $1 <20 and $2==30'|uniq|grep jpg Test.gif

基因组序列

perl -i -pe ‘tr/a-z/A-Z/ unless /^>/’ Hic_genetic.fasta

perl one-liner 参数传递
perl -slane ‘print $p a r a,$ _’ – -para='Print in Perl: ’ in.file

perl的数组内插中进行运算，在简单的语句中还是比较有用
perl -le “while(<*>){print qq{mv $_ @{[s/\s/_/r]}}}”

mv Aanans comosus Aanans_comosus

mv Amborella trichopoda Amborella_trichopoda

当然以上这个，去除当前文件夹中所有文件名中的空格，直接
perl -le “while(<*>){rename $,s/\s//r}”

for file in *;do mv “$file” perl -e '\$name=shift;print $name=~s/[^0-9a-zA-Z.]/_/gr' "$file";done
linux手动更新locate数据库命令 --> updatedb

Perl自定义模块的路径包含问题
http://blog.sina.com.cn/s/blog_a577563c0101353j.html
use FindBin qw( $B in); u se l ib "$ Bin/…/lib";

Use TBtools to facilitate the utilization of NGS data.

1. Extract Fasta Record…

2. Conduct gene ontology enrichment

3. Remote Blast and Visualize the result

4. Select Rows from a huge table, like those that contain gene annotations

5. Compare ID set by Venn Graph, 2 set? 3set? Use TBtools

6. Visualize the expression pattern of several gene in a short time

…

Use TBtools !!!
http://cj-chen.github.io/tbtools/TBtools-seminar/

快速批量提取fasta记录

快速进行本地化blast

2分钟30秒从转录组中鉴定提取感兴趣基因的同源基因（图形化展示比对结果）

sort -k8,8 -k9n,9
http://www.cnblogs.com/longjshz/p/5797933.html

小功能的脚本写了就丢，
不如复制命令行方便，
随取随用 – CJ chen 20140717
（注：window下，one-liner应换两边单引号为双引号,同时，少数命令暂仅限于linux下使用）—随时接收修改意见

最新命令行增加时间： 2021年07月11日17:08
最近命令行优化时间：2015年03月28日11:05
最近BUG修复时间：2016年04月12日09:20

使用tips!!!
1.使用ctrl+F，复制标题，直接检索日志
2.批量处理文件时在linux下使用以下命令
for batch in *.suffix;do 本日志中的命令行 > ${batch}_suffix2 &; done
#解释，由于使用了perl命令行，单引号避免了shell变量内插，
将任务挂至后台，可批量同时运行，实现批处理，去掉该符号，则可逐个文件处理，
使用了共同后缀名，文件名通配
亲测同时清理20个fq文件，哈哈

#目录 ----tip：

-1.查看fastq文件读段平均读长、最大读长、最短读长
0.perl命令行粗暴多文件并行处理（每个线程处理一个文件）
1.从fasta文件中提取特定的某个序列(记录) —好用界面版请见空间日志《Grape_one_fasta.jar》
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输入文件
63. 提取任意组装结果最长转录本(so-called Unigenes)或者CDS预测结果中最长序列
64. 表格类数据，以某一列为keys组成的Group中仅保留其对应某属性（另一列）中值最大的一类
65. 小文件行随机化
66. 打印匹配行及其前’指定数目’行
67. 打印匹配行及其后’指定数目’行
68. -n的多个文件区别对待
69. 按照列名提取文件多列
70. 批量提取多个序列多个区段
71. 输出fasta文件每个序列对应的长度 ID\tLength\n
72. 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比对结果
104. 统计fastq文件的碱基数以及Q20和Q30的碱基数
105. 基于NCBI的GI列表或者Accession列表获取相应物种信息
106. 依据两个tab分割表格的第一列合并两个表格
107. 字符串去冗余（压缩相邻所有相同字符为1个，并输出个数）
108. 整理intron大小，从gff3文件
109. 等分Fasta文件为指定个数
110. 等分表格文件，按行等分为指定个数文件
111. Fastq文件，对reads进行trimming
113. 格式话Fasta序列为序列每行60个碱基
114. Fasta ID简化及信息提取
115. 整理cd-hit-est结果为 gene2trans map 用于RSEM
116. 输出 00001 到 10000
117. 表格转置
118. fasta2table
119. 反向互补子函数
120. 查看bam文件中Deletion分布
121. 对矩阵各列进行求和
122. 对gff文件进行 group // 用于基于位置去重，或者是给 Overlap的record 分组，方便处理
123. 将Sanger测序结果合并，.seq 文件都没有header
124. 解析gff文件的注释信息（MSU Rice，并转换其中的 URI编码）
125. 整理NCBI的Genome database上下载的gff文件中mRNA到protein的ID map
126. 统计fasta文件整体的碱基分布 ATCG…
127. 解析KEGG Blast2KO 的文本 - KEGG htext 格式得到注释表格
128. 过滤和筛选WGCNA输出网络表格的模块信息
129. 使用ID映射列表，替换目标文件中所有目标文本
130. 简化Fasta的ID，去除ID行空格之后的无用信息
131. 批量重命名目录下的文件（由于文件名的部分匹配和替换-基于一个列表）
132. Gff3文件中提取ID映射关系
133. 解析cd-hits-est 的聚类结果为tab (方便后续分析)
134. 整理PlantCARE的预测结果
135. 纯粹perl批量重命名 (适用所有系统，rename)
136. 对行进行Bootstrap抽取（即指定输出行数的可重复抽取）
137. 从eggNOGmapper注释结果中提取GO注释信息

#############################

查看fastq文件读段平均读长、最大读长、最短读长

##############################
#查看全读段数据的平均读长，最大读长，最短读长 ===用以估算数据量大小，查看是否可能是原始数据
perl -ne ‘BEGIN{ $min = 1 e 10;$ max=0;}next if ( $.$ read_count++; $cur_length=length($ _); $total_length+=$ cur_length; $min =$ min> $cur_length?$ cur_length: $min;$ max= $ma x <$ cur_length? $cur_length:$ max;END{print qq{Totally $read_count reads\nMAX length is $max bp\nMIN length is $KaTeX parse error: Undefined control sequence: \nMean at position 8: min bp \̲n̲M̲e̲a̲n̲ ̲length is },$ total_length/$read_count,qq{ bp\n}}’ plum_T1_F.fq

##########################不要让文件数超过线程数，囧

perl命令行粗暴多文件并行处理（每个线程处理一个文件）

##########################
perl -e ‘while(<*.sra>){system (“./tools/fastq-dump $_ &”)}’

####################################

从fasta文件中提取特定的某个序列(记录)

####################################
#注意：此命令使用了正则，但因仅仅提取一个记录，与hash差别不大

perl -e '$file=shift;$record=shift;$record=quotemeta($record);open SEQ,$file or die qq/can_t open file\n/;while(<SEQ>){next unless /$record\b/;print;while(<SEQ>){exit if />/;print}}' <fasta_file> [id]
# 更简单的命令为
perl -0076 -ane 'chomp;print qq{>$_} if $F[0] eq qq{gene_id}' out.fa
###############################
#  从fasta文件中批量提取序列(记录) 
###############################
# 最新版... 更短了
# 批量提取fasta
perl -lne 'if($switch){if(/^>/){$flag=0;m/^>?(\S+).*?$/;$flag=1 if $need{$1};}print if $flag}else{m/^>?(\S+).*?$/;$need{$1}++}$switch=1 if eof(ARGV)' Clean.ids Unigenes.fasta
# 批量过滤fasta
perl -lne 'if($switch){if(/^>/){$flag=1;m/^>?(\S+).*?$/;$flag=0 if $need{$1};}print if $flag}else{m/^>?(\S+).*?$/;$need{$1}++}$switch=1 if eof(ARGV)' NoPlant.ids Unigenes.fasta 
# 旧版本 太长了
perl -e '$id=shift;open ID,qq{<},$id or die qq{Can_t read id file};while(<ID>){chomp;($id_need=$_)=~s/^>?(\S+).*/$1/;$need{$id_need}++;}$fasta=shift;open FA,qq{<},$fasta or die qq{Can_t read fa file};while(<FA>){chomp;if(/^>/){$flag=0;($id_found=$_)=~s/^>(\S+).*/$1/;if ($need{$id_found}){print $_,qq{\n};$flag=1;}}else{if ($flag){print $_,qq{\n};}}}' id_per_line file.fa

# 对于部分文件 id 隐藏得比较深，此时，要用强大的正则匹配... 但速度必然比上面hash的慢，且列表越大越慢
 perl -e "$id=shift;open ID,qq{<},$id or die qq{Can_t read id file};while(<ID>){chomp;push @need,$_;}$fasta=shift;open FA,qq{<},$fasta or die qq{Can_t read fa file};while(<FA>){if(/^>/){$flag=0;$id_found=$_;if (grep {$id_found=~/\b$_\b/} @need){print $_;$flag=1;}}else{if ($flag){print $_;}}}" id_per_line file.fa 

 ###########
#  格式转换  #
###########
#Fastq格式转换为fasta格式http://
perl -ne '$linenum=$.%4;if($linenum==1){s/^@/>/;print;}if($linenum==2){print}'  <fastq.file> > <fasta.file_name>

######################################
#  常规fasta文件去格式为一行id一行seq
#####################################
#  20150105优化，正则匹配 $seq=~s/\r?\n(.)/$1/g其中的(.) perfect ! 一行搞定
perl -ne 'if (/^>/){$seq=~s/\r?\n(.)/$1/g;print $seq;$seq=q{};print;}else{$seq.=$_}END{$seq=~s/\r?\n(.)/$1/g;print $seq;}' in.fa > out.fa

###################################
#    快速批量提取读段文件的指定序列  (也可用于去格式的fasta文件)
#####################################
#从极其简化的fa文件(一行id，一行seq中快速批量提取序列)----优化，将文件的id中无谓的部分暴力去除，方便匹配
perl -e '$reads=shift;$list=shift;open LIST,$list or die "can_t find file";while(<LIST>){chomp;s/^>//;s/\s+$//;$need_id{$_}=1;}open READ,$reads or die "Can_t open reads";while(!eof(READ)){$id=<READ>;$seq=<READ>;chomp $id;$id=~s/^>//;$id=~s/(\S+?\b).*/$1/;print ">",$id,"\n",$seq if ($need_id{$id})}' <fasta.file> <id.list>
#注意到这一命令行仅仅用于id一行，序列一行的fa格式，本来就是为了debug用的


##############
#  读段个数统计 #
##############
#fatsq文件读段数目统计
wc -l <fastq.file> |perl -ane 'print "$F[0]/4\n"'|bc
   也可以直接 grep '@' <fastq.file> | wc -l
 
#fasta文件序列个数统计
 grep '>' <fasta.file> | wc -l

######################################
# fastq质量值格式转换 Phred+64转Phred+33  #
######################################
perl -e 'while(<>){$line=$_;unless($.%4){$line="";chomp;@asc=(split(//,$_));$line.=chr(ord($_)-31) foreach @asc;$line.="\n"}print $line}' <fastqfile phred+64> > <outfile_name>
# 改进版
... 懒得写

 
##################
#  fastq 5'端trimming  #
##################
 #规定的长度短于设定值，则不切，直接修改$length=? 
#切除的碱基数目，直接修改$trimmed_length=?
perl -e 'BEGIN{$length=20;$trimmed_length=8}while(<>){$cont=$_;unless($.%2){unless(length($cont)<$length){$cont=substr($cont,$trimmed_length)}}print $cont}' <fastq_file> > <fastq_file_trimmed>

#########################################
#    去除低质量值碱基数量高于N个的reads---用于phred+33数据
########################################
#多强大的clean工具，都不能适用于一切数据，一定要手工检查，确保工具盲点，和参数设置错误
#从质量值格式为phred+32的fastq文件中，去除所有含有质量值低于20(包括20)超过5个(包括5个)的读段
perl -e '$reads=shift;open READ,$reads or die "Can_t open reads";while(!eof(READ)){$id=<READ>;$seq=<READ>;$add_label=<READ>;$qual=<READ>;chomp($qual);@asc=(split(//,$qual));$low_count=0;foreach(@asc){$asc_number=ord($_)-32;$low_count++ if ($asc_number<=20)}print $id,$seq,$add_label,$qual,"\n" unless ($low_count>=5)}' <fastq.file> > <result.file>

###################
#除读段序列含未知碱基N超过一定比例的读段
######################
#去除N含量多余10%的读段
perl -e 'BEGIN{$drop_N_percent=0.10}$reads=shift;open READ,$reads or die "Can_t open reads";while(!eof(READ)){$id=<READ>;$seq=<READ>;chomp $seq;$add_label=<READ>;$qual=<READ>;$N_count=($seq=~tr/N/N/);print $id,$seq,"\n",$add_label,$qual unless ($N_count/length($seq)>=$drop_N_percent)}' <fastq.file> > <result.file>


############################20140807 bug 已修复，原版会切掉每个读段最后一个碱基
#     切除读段两端质量值低于给定阈值的部分并丢弃长度低于给定值的记录
#################################20140831 修正 并 增加双端版本 
#  20140921之前控制长度少了一个碱基，bug已修复
# 20150809 之前碱基数值算错啦，写成了-32 实际上，应该是 -33
#修改低质量值碱基质量值阈值$cut_qual=20，修改长度阈值$cut_length, 低于此长度的读段将被丢弃
perl -e 'BEGIN{$cut_qual=20;$cut_length=25}$reads=shift;open READ,$reads or die "Can_t open reads";while(!eof(READ)){$id=<READ>;$seq=<READ>;chomp $seq;$add_label=<READ>;$qual=<READ>;chomp $qual;@asc=(split(//,$qual));for($i=0;$i<=$#asc;$i++){last unless ord($asc[$i])-33<$cut_qual};next if ($i==$#asc);for($j=$#asc;$j>=0;$j--){last unless ord($asc[$j])-32<$cut_qual};next unless ($j-$i+1>=$cut_length);if($j==$#asc){print $id,substr($seq,$i),qq/\n/,$add_label,substr($qual,$i),"\n";}else{print $id,substr($seq,$i,$j-$i+1),qq/\n/,$add_label,substr($qual,$i,$j-$i+1),"\n";}}'

#20140831，有两同学需用双端版本，应保证正反向对齐（如原始下机数据） BUG已修，20140921 BUG修复 20150809
perl -e 'BEGIN{$cut_qual=20;$cut_length=5}$reads1=shift;open READ1,$reads1 or die "Can_t open reads";$reads2=shift;open READ2,$reads2 or die "Can_t open reads";$out1=shift;$out2=shift;open OUT1,q{>},$out1 or die "Can_t output reads";open OUT2,q{>},$out2 or die "Can_t output reads";while(!eof(READ1)){$id1=<READ1>;$seq1=<READ1>;chomp $seq1;$add_label1=<READ1>;$qual1=<READ1>;chomp($qual1);@asc1=(split(//,$qual1));$id2=<READ2>;$seq2=<READ2>;chomp $seq2;$add_label2=<READ2>;$qual2=<READ2>;chomp($qual2);@asc2=(split(//,$qual2));for($i1=0;$i1<=$#asc1;$i1++){last unless ord($asc1[$i1])-33<$cut_qual};next if ($i1==$#asc1);for($i2=0;$i2<=$#asc2;$i2++){last unless ord($asc2[$i2])-33<$cut_qual};next if ($i2==$#asc2);for($j1=$#asc1;$j1>=0;$j1--){last unless ord($asc1[$j1])-33<$cut_qual};next unless ($j1-$i1+1>=$cut_length);for($j2=$#asc2;$j2>=0;$j2--){last unless ord($asc2[$j2])-33<$cut_qual};next unless ($j2-$i2+1>=$cut_length);if($j1==$#asc1){print OUT1 $id1,substr($seq1,$i1),qq/\n/,$add_label1,substr($qual1,$i1),"\n";}else{print OUT1 $id1,substr($seq1,$i1,$j1-$i1+1),qq/\n/,$add_label1,substr($qual1,$i1,$j1-$i1+1),"\n";}if($j2==$#asc2){print OUT2 $id2,substr($seq2,$i2),qq/\n/,$add_label2,substr($qual2,$i2),"\n";}else{print OUT2 $id2,substr($seq2,$i2,$j2-$i2+1),qq/\n/,$add_label2,substr($qual2,$i2,$j2-$i2+1),"\n";}}' forward.fq reverse.fq out_forward.fq out_reverse.fq

################################
#    去除低质量值碱基(Q<给定值)所在比例高于(P大于给定值)的读段---用于phred+33数据
###############################
#直接小修上一个命令行
#直接修改$low_qual=20的值指定低质量碱基阈值，修改$drop_percent=0.10,指定低质量碱基最高所占比例阈值（此处10%）
perl -e 'BEGIN{$low_qual=20;$drop_percent=0.10}$reads=shift;open READ,$reads or die "Can_t open reads";while(!eof(READ)){$id=<READ>;$seq=<READ>;$add_label=<READ>;$qual=<READ>;chomp($qual);@asc=(split(//,$qual));$low_count=0;foreach(@asc){$asc_number=ord($_)-33;$low_count++ if ($asc_number<=$low_qual)}print $id,$seq,$add_label,$qual,"\n" unless ($low_count/@asc>=$drop_percent)}' <fastq.file> > <result.file>


##################
# DNA序列转mRNA序列
##################
#DNA->mRNA
perl -npe 'unless(/^>/){tr/ATCG/TAGC/;s/T/U/g}' <fasta.file>

###########################################
#  perl脚本windows和linux间切换，导致程序运行终止 #
###########################################
报错信息为：No such file.....
解决办法，备份同时去除win下面回车符
 perl -i -pe 's/\r//' getFileNames.pl 
 
 #########################
#  window下打印前10行
#########################
perl -pe "exit if $.>10" in.file
# 打印最后10行
perl -ne "push @a,$_;shift @a if (@a>10)};print @a;{" in.file
 
######################
#  生成批处理用file_list
######################
生成没有后缀的文件名列表-----修改其中quotemeta(".fastq"),即可
perl -e 'BEGIN{$old_suffix=quotemeta(".fastq");}while(<*$old_suffix>){s/$old_suffix//;print $_,"\n"}' > your_list


######################
#    fastq中提取特征读段
###########################
#以下例子是提取可能的polyA读段。。。缺点在于，根本区分不了插入片段过小的读段
 perl -e '$reads=shift;open READ,$reads or die "Can_t open reads";while(!eof(READ)){$id=<READ>;$seq=<READ>;$add_label=<READ>;$qual=<READ>;print $id,$seq,$add_label,$qual if ($seq=~/A{18}$/)}' <fastq.file> > <result.file>

#GC含量在某一区间的读段
perl -e '$reads=shift;open READ,$reads or die "Can_t open reads";while(!eof(READ)){$id=<READ>;$seq=<READ>;$add_label=<READ>;$qual=<READ>;chomp $seq;$GC=$seq=~tr/GC/GC/;print $id,$seq,qq{\n},$add_label,     $qual if ($GC>=55&&$GC<=63)}' white_F_final_highqual_clean.fq > white_GC_55_63.fq

# 小RNA数据过滤长度  只保留长度在35bp到53bp的读段
perl -e '$reads=shift;open READ,$reads or die "Can_t open reads";while(!eof(READ)){$id=<READ>;$seq=<READ>;$add_label=<READ>;$qual=<READ>;chomp $seq;print $id,$seq,qq{\n},$add_label, $qual if (length($seq)>=35&&length($seq)<=53)}' white_F_final_highqual_clean.fq > white_length_35_53.fq

#    提取指定id（除了主要barcodes外的fq记录） 20150403 将unless修改为if，之前写错了
perl -e '$reads=shift;open READ,$reads or die "Can_t open reads";while(!eof(READ)){$id=<READ>;$seq=<READ>;$add_label=<READ>;$qual=<READ>;print $id,$seq,$add_label,$qual if ($id=~/ACAGTG$/)}'

#############################
#    fasta格式CDS转为aa（必须有终止密码子）
##################################
perl -MBio::SeqIO -MBio::Tools::CodonTable -e '$file=shift;$fasta=Bio::SeqIOnew(-file=>$file);while($seq=$fasta->next_seq){print $seq->id,"\n",$seq->translate(-complete=>1)->seq,"\n"}' <fasta.file>
#必须有终止密码子 

############################
#    window下面模拟cut命令-提取文本第二列
##########################
#提取第二列
perl -ane "print $F[1],qq/\n/"

##############################
#    window下合并多个fasta文件
###########################
perl -e "while(<*.fa>){system qq{type $_ >> combined.fa}}" 

######################################
#     window下提取匹配到某一模体的fasta序列
#####################################
# window下常规fasta文件去格式为一行id一行seq
perl -pe "s/\r?\n//;if(/^>/){s/^/\n/;s/$/\n/}END{print qq/\n/}" <fasta.file> > <result.fasta.file>
perl -ne "print unless /^\s+$/" <result.fasta.file> > <final.result.fasta.file>
# 提取匹配到某一模体的序列
perl -e "$file=shift;$pattern=shift;chomp $pattern;$pattern=quotemeta($pattern);open IN,$file or die qq/can_t open the sequence file/;while(!eof(IN)){$id=<IN>;$seq=<IN>;print if ($seq=~/$pattern/)}" <final.result.fasta.file> YOURPATTREN > <result_file_name>

# 一行版本
perl -e 'my $label_pattern=shift;my $file=shift;my $NeedSeqPattern=shift;my $ref=&Single_Label_Reader($file,$label_pattern);my $content;while($content=&$ref){my ($label,@content)=@$content;chomp ($label,@content);my $seq=join q{},@content;print qq{$label\n$seq\n} if ($seq=~/$NeedSeqPattern/i)}sub Single_Label_Reader {my $file=shift;my $label_pattern=shift;die "Error:$file unreachable.\n" unless (-s $file);open IN,q{<},$file or die "Error:$file unreadable.\n";my $temp1=<IN>;my $temp2=$temp1;my @temp;return sub {return q{} if (eof(IN));while(!eof(IN)){$temp1=$temp2;@temp=();while(<IN>){if (/$label_pattern/){$temp2=$_;last;}push @temp,$_;}return [$temp1,@temp];}}}' '^>' in.fasta ATATATATA


##################################
#    提取人类基因组注释文件rRNA注释
##################################
#有朋友使用cuf分析人类转录组数据的时候，需要使用-M选项，掩盖掉rRNA信息
perl -ne 'print if /rRNA/' <hg75.gtf>
#但这样做也会把其他rRNA,如mt上的rRNA注释提取出来
使用
perl -ne 'print if /\brRNA\b/' <hg75.gtf>
即可避免

#####################################
#    对sort | uniq -c | 的结果频次由高到低排序，有大用
######################################## 
#生成结果直接用于blast   #bug已修复 #bug已修复 20140909
grep -P 'AGCCTAC' sample_F.nolowbase.clean|sort|uniq -c|perl -ne '($count,$seq)=split;$count{$seq}=$count;END{foreach (sort {$count{$b}<=>$count{$a}} keys %count){print $_,"\t",$count{$_},"\n"}}'|head|cut -f1|perl -ne 'print ">$.\n$_"'


############################
#    fasta格式的DNA序列反向互补---仅使用于一行id一行序列的fa文件
########################
 perl -ne 'if (/^>/){print}else{s/\r?\n//;tr/atcgATCG/tagcTAGC/;@a=split(q//,$_);print reverse @a;print qq/\n/}'


######################
#    一行id一行序列的fa文件格式化为一行id多行序列
#############################
perl -ne 'BEGIN{$length=60;}if(/^>/){print}else{$sequence=$_;for($pos=0;$pos<length($sequence);$pos+=$length){print substr($sequence,$pos,$length),"\n"}}' Ananas.Unigene.filtered.fa|perl -ne 'print unless /^\s+$/' > Ananas.Unigene.filtered.fmt.fa

###################---使用Bowtie对读段进行深度清理后，unmapped.fastq顺序乱，需要整理，对应双端
#     按fastq文件标签名对读段顺序进行排序---待优化版(目前需10倍于文件大小的内存)
########################-----具体优化思路，使用Tie::File模块，每个读段仅保存其id即id对应的行号
perl -e '$reads=shift;open READ,$reads or die "Can_t open reads";while(!eof(READ)){$id=<READ>;$seq=<READ>;$add_label=<READ>;$qual=<READ>;$info{$id}=[$seq,$add_label,$qual]}for(sort keys %info){print $_,@{$info{$_}}}'

#####################---注意-对于双端文件，应在原始数据或者已配对数据上同时修改，否则后期无法对齐
#替换fq或fa文件记录的id为指定形式
################################
#修改$newid=qq/my_id/即可用于修改fq文件的id,修改$.%4==1为$.%2==1即可用于修改fa文件的id
perl -ne 'BEGIN{$newid=qq/my_id/}if($.%2==1){print $newid,qq/-/,++$i,qq/\n/}else{print}'

#Trinity组装结果，去除剪接本信息，仅保留至基因水平的id信息
perl -pe 's/(>c\d+_g\d+).*/$1/ if /^>/' Longest_trans.fasta

######################
#    提供一个序列名列表逐一替换fasta记录的id    --window下，linux下请将单引号改为双引号
#############################
#...序列名列表，一行一个
perl -e "$fasta=shift;$name=shift;open NAME,$name or die;open FASTA,$fasta or die;while(<NAME>){chomp;$name=$_;while(<FASTA>){if (/^>/){print qq{>$name\n};last};print}}END{while(<FASTA>){print}}" <file.fa> <name.txt>

################################
#   根据NCBI gene id 即gi号获取GeneBank上的序列
##############################
#只有一个gi号,修改其中gi号即可
perl -MBio::DB::GenBank -e '$gb = Bio::DB::GenBank->new();$seq = $gb->get_Seq_by_gi('433680129');print q{>},$seq->accession_number,qq{\n},$seq->seq,qq{\n}'
#批量获取，则提供一个list文件，每行一个gene_id
perl -MBio::DB::GenBank -ne 'chomp;push @my_gene_ids,$_;END{$gb = Bio::DB::GenBank->new();$seq = $gb->get_Stream_by_id(\@my_gene_ids);while(my $GOI=$seq->next_seq){print q{>},$GOI->accession_number,qq{\n},$GOI->seq,qq{\n}}}' <Your_gene_id_list_one_id_per_line> 

#上面的命令不够爽，构建进化树有时候要nucl序列，挑战一下
############################
#根据蛋白gene_id或accession获取其Genebank上的核苷酸序列
###########################
#--不确保绝对正确或可用 #其中关键在于该样本具体链接到几个库
#试验了一下，
#常规序列基本是1个库，其中{'dblink'}->[1]应该修改为{'dblink'}->[0]
#由高通量测序所得数据基本有两个库,其中{'dblink'}->[1]

##以下为获取单个蛋白对应NCBI上核酸序列的命令,直接修改其中ID即可
perl -MBio::DB::GenBank -e '$gb = Bio::DB::GenBank->new();$seq = $gb->get_Seq_by_gi('XP_007097302');$dblink=$seq->annotation->hash_tree->{'dblink'}->[1];$myGirl= $dblink->primary_id;$final_seq=$gb->get_Seq_by_gi($myGirl);print q{>},$final_seq->accession_number,qq/\n/,$final_seq->seq,qq/\n/'

##以下批量获取蛋白对应NCBI上核酸序列的命令
perl -MBio::DB::GenBank -ne 'chomp;push @my_gene_ids,$_;END{$gb = Bio::DB::GenBank->new();@new_id=map{$seq = $gb->get_Seq_by_gi($_);$dblink=$seq->annotation->hash_tree->{'dblink'}->[1];$myGirl= $dblink->primary_id}@my_gene_ids;$final_seq=$gb->get_Stream_by_gi(\@new_id);while(my $GOI=$final_seq->next_seq){print q{>},$GOI->accession_number,qq{\n},$GOI->seq,qq{\n}}}' <列表文件_每行一个id>

############################
#     比较字符串中两个单字符的频次(比如投票0,1或方向F,R)
##########################
 perl -ne 'print tr/R/R/>tr/F/F/?qq{R\n}:qq{F\n}'

################################################
#    从每个记录行数确定的文件中随机抽取不重复的N个记录
#################################################
#提供三个参数，
#1.文件，
#2.每个记录所占行数（必须每个一样，比如fastq文件 参数为4 ），如果要提取fa文件，可使用上文flat fasta的命令，再使用此命令
#3.所需记录的总数，比如我要1000个读段的数据，那么就给 1000
#仅适用于linux下，因需要文件具体行数
perl -e '$file=shift;$lines_per_record=shift;$sample_num=shift;$total_line=`wc -l $file`;$total_line=~s/^(\d+).*/$1/;chomp $total_line;$total_record=$total_line/$lines_per_record;for(1..$sample_num){$int=int(rand($total_record)+1);unless ($hash{$int}){$hash{$int}=1}else{redo}};open FILE,q{<},$file or die qq{Can_t open $file};while(!eof(FILE)){for(1..$lines_per_record){$temp_cont=<FILE>;push @temp_record,$temp_cont;}$count++;if($hash{$count}){print @temp_record;}@temp_record=();}' <your_file> [number_of_lines_per_record] [number_of_records_tobe_extracted]


#################
#    有同学想知道比对上的读段在genome上正反链的分布情况
####################    
perl -ane 'if($F[1]&0x0010){print "$F[0]\t-\n"}else{print "$F[0]\t+\n"}' your_sam_file


######################
#    去除全读段所有碱基质量值均低于某个阈值（如20）的读段(支持单端和双端数据)
############################
#不用考虑，这些读段绝对不能要
#对于单端数据
perl -e 'BEGIN{$low_qual=20;}$reads=shift;open READ,$reads or die "Can_t open reads";while(!eof(READ)){$id=<READ>;$seq=<READ>;$add_label=<READ>;$qual=<READ>;chomp($qual);@asc=(split(//,$qual));$don_t_throw=0;foreach(@asc){$asc_number=ord($_)-32;if($asc_number>=$low_qual){$don_t_throw++;last}}print $id,$seq,$add_label,$qual,"\n" if ($don_t_throw)}' <fastq.file> > <result.file>

#对于双端数据--注意，利用了双端数据必定配对存在
 perl -e 'BEGIN{$low_qual=20;}$read1=shift;$read2=shift;$out1=shift;$out2=shift;open READ1,$read1 or die "Can_t open reads1";open READ2,$read2 or die "Can_t open reads2";open OUT1,q{>},$out1 or die "Can_t output reads1";open OUT2,q{>},$out2 or die "Can_t output reads2";while(!eof(READ1)){$id1=<READ1>;$seq1=<READ1>;$add_label1=<READ1>;$qual1=<READ1>;chomp($qual1);$id2=<READ2>;$seq2=<READ2>;$add_label2=<READ2>;$qual2=<READ2>;chomp($qual2);$don_t_throw=0;@asc1=(split(//,$qual1));foreach(@asc1){$asc_number1=ord($_)-32;if($asc_number1>=$low_qual){$don_t_throw++;last}}if($don_t_throw){@asc2=(split(//,$qual2));foreach(@asc2){$asc_number2=ord($_)-32;if($asc_number2>=$low_qual){$don_t_throw++;last}}}next unless ($don_t_throw==2);print OUT1 $id1,$seq1,$add_label1,$qual1,"\n";print OUT2 $id2,$seq2,$add_label2,$qual2,"\n"}' <fastq1.file> <fastq2.file> <result1.file> <result2.file>

####################
#    借用pileup文件直接统计测序数据在各染色体上的分布，注意 mpileup -A --尚未判定意义，可能是覆盖率分布
######################
perl -ane '$cover_count{$F[0]}+=$F[3];END{$total+=$_ foreach (values %cover_count);foreach(keys %cover_count){print $_,"\t",$cover_count{$_}/$total,"\n"}}' ab2.pileup

#####################
#    查看sam中uniq mapped比率
########################
cat sampe_a700.sam|perl -ane 'next if /^@/;$count++;;$uniqmap{$F[11]}++;END{for(keys %uniqmap){print $_,"\t",$uniqmap{$_},"\t",$uniqmap{$_}/$count,"\n"}}' >XT_A_stat_of_sampe_a700.sam

###########
#    查看sam中编辑距离分布
#################
 cat sampe_a700.sam | perl -ane 'next if /^@/;$count++;;$uniqmap{$F[12]}++;END{for(keys %uniqmap){print $_,"\t",$uniqmap{$_},"\t",$uniqmap{$_}/$count,"\n"}}' > NM_i_stat_of_sampe_a700.sam

##############
#  统计各行平均值或各列平均值
############
#统计各行平均值
perl -ane '$sum+=$_ for(@F);print $sum/@F,"\n"'
#统计各列平均值
perl -ane 'for(0..$#F){$sum[$_]+=$F[$_]};END{print $_/$.,"\n"for(@sum)}'

#统计各行的平均值，如果第一列是 基因ID的话 
perl -F'\t' -lane 'print $F[0],qq{\t},eval((join qq{+},@F[1..$#F]))/$#F'


################
#    将fa文件(尤其基因组文件)分成每个记录一个文件(要求一行id一行seq,见25)
####################
perl -e 'while(<>){$label=$_;$seq=<>;$file=$label;$file=~s/^>//;chomp($file);open OUT,q{>},$file;print OUT $label,$seq}' yourchromosome.fa 
#实际上，这类文件，可直接用linux命令，split -l 2 文件    即可，用命令的好处仅仅是win下通用以及分出的文件文件名即id


#############
#    批量重命名 win下修改单引号为双引号
##############
perl -e 'while(<sampe_a700*>){$newname=$_;$newname=~s/sampe_a700/DY1/;rename $_,$newname}'

#############
#    win下批量去除文件夹内所有文件中的数字  linux下修改双引号为单引号
#############
 perl -e "while(<*>){open IN,$_;open OUT,q{>},q{new}.$_;while(<IN>){s/\d//g;print OUT;}}" 


##############
#    统计SAM文件某一标签(BWA结果)
###############
perl -ane '@x1=map {$_ =~ /X1:i:(\S+)/} @F;print $x1[0],"\n" if (@x1);' DY.pre.norm.sam | sort |uniq -c|sort -n > DY.pre.norm.sam.x1_stat

##################
#   提取长度大于1000bp的fa记录
##################
# 更新了
perl -ne 'BEGIN{$id=q{};$seq=q{}}chomp;if(/^>/){print $id,qq{\n},$seq,qq{\n} if (length($seq=~s/\s|-//gr)>=474);$id=$_;$seq=q{}}else{$seq.=$_;}END{print $id,qq{\n},$seq,qq{\n} if (length($seq=~s/\s|-//gr)>=474)}' PMMV.fas|grep -c '^>'
# 简化版 
perl -0076 -ne 'chomp;($id,@seq)=split /\r?\n/,$_,2;print qq{>$id\n@seq} if length(join q{},map {s/\s+|-//gr;}@seq)>=474' PMMV.fas|grep -c '^>
 
#window下请自行修改两端单引号为双引号
perl -ne 'BEGIN{$id=q{};$seq=q{}}chomp;if(/^>/){print $id,qq{\n},$seq,qq{\n} if (length($seq)>=1000);$id=$_;$seq=q{}}else{$seq.=$_;}END{print $id,qq{\n},$seq,qq{\n} if (length($seq)>=1000)}' file.fasta

####################
#    批量提取匹配行（正则匹配，强大）
#####################
#修改后可用与批量提取展平的fasta文件，非常强大的grep  正则匹配是^ \b 非常重要！！！
perl -e '$list=shift;$fpkm=shift;open LIST,$list;@ids=<LIST>;chomp @ids;open FPKM,$fpkm;while(<FPKM>){$line=$_;print $line if (map {$line=~/^$_\b/} @ids)}'  complete.orfs.ids PINEAPPLE_RSEM.isoforms.results > complete.orfs.fpkm


#################
#  fasta中有相同id，增加后缀方便blast建库
####################
perl -lpe 'next unless /^>/;$hash{$_}++;if($hash{$_}>1){$_.=qq{_$hash{$_}}}'


##################
#    多个列表文件，比如gene_ids，取样品特异gene_id 
###############
perl -e 'for(1..3){push @files,shift}for(@files){$sample=$_;open IN,$_;while(<IN>){$counts{$_}++;push @{$samples{$sample}},$_}}for(keys %samples){$sample=$_;for(@{$samples{$sample}}){print "$sample\t$_" if ($counts{$_}==1)}}' *.Genes > single_genes.xls


###############
#  直接统计一个序列的GC含量
###############
perl -e '$seq=shift;chomp $seq;$GC=$seq=~tr/GC/GC/;print $GC/length($seq),qq{\n}' CTTTCATCTTCCCTGCAAAGAGTGTCAATAGATCTCCATGGCCGCGATTTCT

################
#    直接连接几个序列并将小写转换成大写
#####################
perl -e '$seq=join ("",@ARGV);$seq=~tr/atcg/ATCG/;print $seq,qq{\n}' cttggaaa tgcagcactg tctggtacgt

##################
#    序列贪吃蛇
###################
#撰写原因：
# 下机数据有overrepresent序列或者子串，有可能是污染读段或者接头等，使用本命令可将其尽可能全面地抓取出来，并用于blast检测
#参数解释：
# $index_length为每次检索序列的长度，默认为25
# $seek_time为遍历文件次数，一般遍历3次，有时候我会遍历5次
# fastq_file 为提供输入的fastq文件
# seq_pattern 为提供的用于一开始进行检索的字串
#结果解释：
# 返回一个由提供子串开始拼接而成的序列，可直接用于blast
perl -e 'BEGIN{$index_length=25;$seek_times=3;}$file=shift;$seq=shift;$total_seq.=$seq;open FQ,$file or die qq{Please input a FASTQ file};while(!eof(FQ)){while(<FQ>){next unless ($.%4==2);if (/$seq/){$line=$_;last;};}if(eof(FQ)&&$seek_times){$seek_times--;print qq{SEEK_TIMEs_REMAIN:},$seek_times,qq{\n};seek FQ,0,0;}if ($line){chomp $line;$pos=index($line,$seq);$new_part=substr($line,$pos+length($seq));print qq{FOUND:$line\n};print qq{NEWPART:$new_part\n};$total_seq.=$new_part;$seq=substr($line,length($line)-$index_length);print qq{INDEX:$seq\n};}}print qq{TOTAL SEQ:\n$total_seq\n};print qq{You may modify the tail of the total_seq to look for possible mistake made by the perl script and get longer total_seq\n};' fastq_file seq_pattern

#以上为向子串下游拼接，以下版本是向上游拼接
perl -e 'BEGIN{$index_length=25;$seek_times=3;}$file=shift;$seq=shift;$total_seq.=$seq;open FQ,$file or die qq{Please input a FASTQ file};while(!eof(FQ)){while(<FQ>){next unless ($.%4==2);if (/$seq/){$line=$_;last;};}if(eof(FQ)&&$seek_times){$seek_times--;print qq{SEEK_TIMEs_REMAIN:},$seek_times,qq{\n};seek FQ,0,0;}if ($line){chomp $line;$pos=index($line,$seq);$new_part=substr($line,0,$pos);print qq{FOUND:$line\n};print qq{NEWPART:$new_part\n};$total_seq=$new_part.$total_seq;$seq=substr($line,0,$index_length);print qq{INDEX:$seq\n};}}print qq{TOTAL SEQ:\n$total_seq\n};print qq{You may modify the tail of the total_seq to look for possible mistake made by the perl script and get longer total_seq\n};' fastq_file seq_pattern


###########
#    随机提取一定比例的fasta或者fastq记录
############
#
perl -e '$file=shift;$lines_per_record=shift;$sample_percent=shift;$total_line=`wc -l $file`;$total_line=~s/^(\d+).*/$1/;chomp $total_line;$total_record=$total_line/$lines_per_record;for(1..int($sample_percent*$total_line/$lines_per_record)){$int=int(rand($total_record)+1);unless ($hash{$int}){$hash{$int}=1}else{redo}};open FILE,q{<},$file or die qq{Can_t open $file};while(!eof(FILE)){for(1..$lines_per_record){$temp_cont=<FILE>;push @temp_record,$temp_cont;}$count++;if($hash{$count}){print @temp_record;}@temp_record=();}' <your_file> [number_of_lines_per_record] [percentage_to_be_extract:0.0-1.0]

###########
#    单行记录随机分组
###############
#perl 实现单行记录 随机分组
perl -MList::Util=shuffle -e "$per=3;@all=shuffle <>;chomp @all;$end3=int(@all/$per)-1;while(@all){for(0..$end3){push @{$list[$_]},shift @all;}}print join qq{\t},@{$_},qq{\n} for(@list)" 

#############
#    按照fasta序列长度排序 (若是直接使用的同学，请直接重命名fasta文件为fasta，直接运行一下命令即可)
#############
#调试用，专门针对fasta记录的字节记录
perl -e '$file=shift;open IN,$file;$pre=tell(IN);while(<IN>){if(/^>/){print qq{$length\t$pre\n$pre\t};$length=0}else{chomp;$length+=length($_)};$pre=tell(IN);}END{print qq{$length\t$pre\n}}' fasta|sed '1d' > pos
----------------------
0       12      21
21      6       37
37      9       55
----------------------
#输出fasta序列长度排序好的pos
perl -e '$file=shift;open IN,$file;$pre=tell(IN);while(<IN>){if(/^>/){print qq{$length\t$pre\n$pre\t};$length=0}else{chomp;$length+=length($_)};$pre=tell(IN);}END{print qq{$length\t$pre\n}}' fasta|sed '1d'|sort -k2 -n|perl -ne '@F=split /\s+/;print qq{$F[0]\t$F[2]\n}' > pos
#使用之前的命令，依据起始字节位置和终止字节位置直接操纵文件，提取序列
perl -e '$file=shift;$pos=shift;open IN,$file;open POS,$pos;while(<POS>){chomp;($start,$end)=split /\s+/,$_;seek(IN,$start,0);read(IN,$content,$end-$start);print $content}' fasta pos >final_fasta

###########
#    双标签区段提取 (可以很方便的用于去除不需要的区块)
###########
# 下一行立即检验右操作数的玩法 ..
perl -ne 'print if /Hsp_hit-from/../Hsp_hit-to/;' unigene.FASTA_blastnr_result.xml
# 还有另外一种玩法，下一行不立即检验右操作数目的 ...
perl -ne 'print if /<\/?Iteration>/.../<\/?Iteration>/' unigene.FASTA_blastnr_result.xml

#比如去除blast 前面的header
perl -ne 'next if /xml version/../BlastOutput_iterations/;print' blast.xml > out.xml

#############
#   批量从uniprot上下载序列
###############
# 准备好每个id一行的列表
 # 使用以下命令完成下载
perl -MLWP::Simple -ne 'chomp;$url="http://www.uniprot.org/uniprot/$_.fasta";print get($url);' ids

###############
#    准备trimmomatic所需的adapter.fa文件
###############
# perl -e 命令 正向接头 反向接头
perl -e '$for=shift;$rev=shift;($for_rev=join q{},(reverse split q{},$for))=~tr/ATGC/TACG/;($rev_rev=join q{},(reverse split q{},$rev))=~tr/ATGC/TACG/;print ">PrefixPE/1\n$for\n";print ">PrefixPE/2\n$rev\n";print ">PE1\n$for\n";print ">PE1_rc\n$for_rev\n";print ">PE2\n$rev\n";print ">PE2_rc\n$rev_rev\n"' AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA

###########
#    提取fasta文件特定记录的特定区段
##########
perl -e '$fa=shift;$id=shift;$start=shift;$end=shift;$id=~s/\r?\n$//;open IN,$fa;while(<IN>){next unless /^>/;last if /$id/}while(<IN>){chomp;last if /^>/;$seq.=$_;}print ">$id\n",substr($seq,$start,$end-$start+1),"\n"' NC_001802.fa NC_001802 1089 1203


##############
#  获取GO term Level 2的信息
############### 
## 输出GO的BP MF CC的所有子term 也就是Level 2;
## 对于非模式生物，得到GO注释信息和兴趣基因的id列表之后可以自行经常超几何检验或者卡方检验
## 具体GO注释的 
## 直接百度 molecular function GO
perl -e 'open TEMP,">bp_children.temp";print TEMP "library(GO.db);xx <- as.list(GOBPCHILDREN);xx\$\`GO:0008150\`";close TEMP;$children=`Rscript bp_children.temp`;@GO=map {/(GO:\d+)/g} $children;for(@GO){print "$_\n"}' 2> error 1> BP_level2_ids;perl -e 'open TEMP,">bp_children.temp";print TEMP "library(GO.db);bp_children<-read.table(\"BP_level2_ids\",header=F);bp_children<-as.character(unlist(bp_children));select(GO.db, keys=bp_children, columns=c(\"TERM\",\"ONTOLOGY\"),keytype=\"GOID\")";';Rscript bp_children.temp 2>error 1> BP_level2_info
perl -e 'open TEMP,">bp_children.temp";print TEMP "library(GO.db);xx <- as.list(GOMFCHILDREN);xx\$\`GO:0003674\`";close TEMP;$children=`Rscript bp_children.temp`;@GO=map {/(GO:\d+)/g} $children;for(@GO){print "$_\n"}' 2> error 1> MF_level2_ids;perl -e 'open TEMP,">bp_children.temp";print TEMP "library(GO.db);bp_children<-read.table(\"MF_level2_ids\",header=F);bp_children<-as.character(unlist(bp_children));select(GO.db, keys=bp_children, columns=c(\"TERM\",\"ONTOLOGY\"),keytype=\"GOID\")";';Rscript bp_children.temp 2>error 1> MF_level2_info
perl -e 'open TEMP,">bp_children.temp";print TEMP "library(GO.db);xx <- as.list(GOCCCHILDREN);xx\$\`GO:0005575\`";close TEMP;$children=`Rscript bp_children.temp`;@GO=map {/(GO:\d+)/g} $children;for(@GO){print "$_\n"}' 2> error 1> CC_level2_ids;perl -e 'open TEMP,">bp_children.temp";print TEMP "library(GO.db);bp_children<-read.table(\"CC_level2_ids\",header=F);bp_children<-as.character(unlist(bp_children));select(GO.db, keys=bp_children, columns=c(\"TERM\",\"ONTOLOGY\"),keytype=\"GOID\")";';Rscript bp_children.temp 2>error 1> CC_level2_info

#############
#    单标签语句块读取
#################
# 解析fasta fastq blast结果等等
perl -e 'my $label_pattern=shift;my $file=shift;my $ref=&Single_Label_Reader($file,$label_pattern);my $content;while($content=&$ref){my ($label,@content)=@$content;chomp ($label,@content);......}sub Single_Label_Reader {my $file=shift;my $label_pattern=shift;die "Error:$file unreachable.\n" unless (-s $file);open IN,q{<},$file or die "Error:$file unreadable.\n";my $temp1=<IN>;my $temp2=$temp1;my @temp;return sub {return q{} if (eof(IN));while(!eof(IN)){$temp1=$temp2;@temp=();while(<IN>){if (/$label_pattern/){$temp2=$_;last;}push @temp,$_;}return [$temp1,@temp];}}}' '^>' in.fasta

#  示例，如解析保留一个比对描述和比对结果的pairwise文件，输出有比对结果的query 长度和sub 长度
perl -e 'my $label_pattern=shift;my $file=shift;my $ref=&Single_Label_Reader($file,$label_pattern);my $content;while($content=&$ref){my ($label,@content)=@$content;chomp ($label,@content);my ($query)=map{/Query= (\S+)/}$label;my ($length1,$length2)=map{/Length=(\d+)/g} @content;next unless $length2;my ($subject)=map {/^> (\S+)/} @content;print qq{$query\t$length1\t$subject\t$length2\n}}sub Single_Label_Reader {my $file=shift;my $label_pattern=shift;die "Error:$file unreachable.\n" unless (-s $file);open IN,q{<},$file or die "Error:$file unreadable.\n";my $temp1=<IN>;my $temp2=$temp1;my @temp;return sub {return q{} if (eof(IN));while(!eof(IN)){$temp1=$temp2;@temp=();while(<IN>){if (/$label_pattern/){$temp2=$_;last;}push @temp,$_;}return [$temp1,@temp];}}}' '^Query'  old_100_fasta_2_v2_fasta.pairwise

# 提取任意fa文件中，序列匹配到某模式的所有记录
perl -e 'my $label_pattern=shift;my $file=shift;my $NeedSeqPattern=shift;my $ref=&Single_Label_Reader($file,$label_pattern);my $content;while($content=&$ref){my ($label,@content)=@$content;chomp ($label,@content);my $seq=join q{},@content;print qq{$label\n$seq\n} if ($seq=~/$NeedSeqPattern/i)}sub Single_Label_Reader {my $file=shift;my $label_pattern=shift;die "Error:$file unreachable.\n" unless (-s $file);open IN,q{<},$file or die "Error:$file unreadable.\n";my $temp1=<IN>;my $temp2=$temp1;my @temp;return sub {return q{} if (eof(IN));while(!eof(IN)){$temp1=$temp2;@temp=();while(<IN>){if (/$label_pattern/){$temp2=$_;last;}push @temp,$_;}return [$temp1,@temp];}}}' '^>' in.fasta ATATATATA

##############
#    核酸序列互补配对的子函数
#####################
sub rev_com {my $seq=shift;chomp $seq;my @seq=split qq{},$seq;@seq=reverse @seq;$seq=join q{},@seq;$seq=~tr/ATCG/TAGC/;return $seq;}
perl -e '$dna=shift;print &rev_com($dna),qq{\n};sub rev_com {my $seq=shift;chomp $seq;my @seq=split qq{},$seq;@seq=reverse @seq;$seq=join q{},@seq;$seq=~tr/ATCG/TAGC/;return $seq;}' ATCG
CGAT

#############
#    分隔fa文件 fq文件 genebank文件 为数据小文件
################
#     使用了本文的单标签提取函数，只需要小修，可以将多类文件简单分割
perl -e 'my $label_pattern=q{>};my $file=shift;my $num=shift;for(0..$num-1){open $aa[$_],q{>},qq{$file.$_};}my $ref=&Single_Label_Reader($file,$label_pattern);my $content;while($content=&$ref){my ($label,@content)=@$content;chomp ($label,@content);my $seq=join q{},@content;$fh=$aa[$count++%$num];print $fh qq{$label\n$seq\n};}sub Single_Label_Reader {my $file=shift;my $label_pattern=shift;die "Error:$file unreachable.\n" unless (-s $file);open IN,q{<},$file or die "Error:$file unreadable.\n";my $temp1=<IN>;my $temp2=$temp1;my @temp;return sub {return q{} if (eof(IN));while(!eof(IN)){$temp1=$temp2;@temp=();while(<IN>){if (/$label_pattern/){$temp2=$_;last;}push @temp,$_;}return [$temp1,@temp];}}}' in.fasta 10

############
#    序列格式化成每行等长并打印的子函数
####################
sub print_sequence{my ($sequence,$length)=@_;for (my $pos=0;$pos<length($sequence);$pos+=$length){print substr($sequence,$pos,$length),"\n";}}

##########
# 从公司返还的注释结果中提取query2gi2GO.table -- for blast3go
##############
# 大部分公司返还的是一个后缀名为xlsx或者xls的文件
# 在windows下用excel打开该文件，文件另存为csv格式，逗号分隔
# 在windows下使用以下命令提取(linux下换双引号为单引号)
perl -ne "chomp;$query = (split /,/,$_)[0];print qq{$query};print qq{\t};@GIs=map{/gi\|(\d+)/g}$_;if(@GIs){print join qq{ },@GIs}else{print qq{No_hits}};print qq{\t};@GOs=map {/(GO:\d+)/g}$_;if(@GOs){print qq{[};print join qq{, },@GOs;print qq{]};}else{print qq{No_hits}}print qq{\n}" 
# 手动编辑,查看是否第一行是标题行，去除即可用于blast3go导入

#####
#   blast2go anno文件转换成blast3go输入文件
#######
perl -ne "@F=split /\t/,$_;unless($seen{$F[0]}++){if ($.==1){next}print qq{$tempName\tHITS\t[};print join qq{, },@temp;print qq{]\t\n};@temp=();}$tempName=$F[0];push @temp,$F[1];END{print qq{$tempName\t\t[};print join qq{, },@temp;print qq{]\t\n};}" blastxml.anno > blast3go.currenttable

###########
#   提取任意组装结果最长转录本(so-called Unigenes)或者CDS预测结果中最长序列
##############
# perl 命令 输入文件 提取id标识符的pattern > 输出文件
# 本例使用的数据文件格式为
>comp100107_c0_seq1|m.115111 comp100107_c0_seq1|g.115111 ORF comp100107_c0_seq1|g.115111 comp100107_c0_seq1|m.115111 type:complete len:282 (-) comp100107_c0_seq1:199-1044(-)
....(命令已考虑了多行和单行序列)
>comp100107_c0_seq1|m.115112 comp100107_c0_seq1|g.115112 ORF comp100107_c0_seq1|g.115112 comp100107_c0_seq1|m.1
15112 type:complete len:201 (+) comp100107_c0_seq1:578-1180(+)
....
>......
# 需要自行准备一个pattern ，本例使用的是  '^>(\S+?)\|.*$' 用于捕捉 comp100107_c0_seq1
# 具体命令如下
perl -e 'my $file=shift;my $sigPattern=shift;my $label_pattern=q{^>};my $ref=&Single_Label_Reader($file,$label_pattern);my $content;while($content=&$ref){my ($label,@content)=@$content;$seq=join q{},@content;$seq=~s/[\r\n]//g;;;;($sig=$label)=~s/$sigPattern/$1/;;;;$seenLen{$sig}=length($seq) if $seenLen{$sig}<length($seq);}my $ref=&Single_Label_Reader($file,$label_pattern);my $content;while($content=&$ref){my ($label,@content)=@$content;$seq=join q{},@content;$seq=~s/[\r\n]//g;;;;($sig=$label)=~s/$sigPattern/$1/;;;;print $label,@content if $seenLen{$sig}<=length($seq);}sub Single_Label_Reader {my $file=shift;my $label_pattern=shift;die "Error:$file unreachable.\n" unless (-s $file);open IN,q{<},$file or die "Error:$file unreadable.\n";my $temp1=<IN>;my $temp2=$temp1;my @temp;return sub {return q{} if (eof(IN));while(!eof(IN)){$temp1=$temp2;@temp=();while(<IN>){if (/$label_pattern/){$temp2=$_;last;}push @temp,$_;}return [$temp1,@temp];}}}' CDS_filtered_by_CD_HIT.fasta '^>(\S+?)\|.*$' > CDS_filtered_by_CD_HIT.fasta.relateMaxLen
# 注意，此处提供的模式 不要给　/．．．/　否则出错

# 提取最长转录本，如果最长的有超过1个 仅保留第一个
perl -e 'my $file=shift;my $sigPattern=shift;my $label_pattern=q{^>};my $ref=&Single_Label_Reader($file,$label_pattern);my $content;while($content=&$ref){my ($label,@content)=@$content;$seq=join q{},@content;$seq=~s/[\r\n]//g;($sig=$label)=~s/$sigPattern/$1/;$seenLen{$sig}=length($seq) if $seenLen{$sig}<length($seq);}my $ref=&Single_Label_Reader($file,$label_pattern);my $content;while($content=&$ref){my ($label,@content)=@$content;$seq=join q{},@content;$seq=~s/[\r\n]//g;($sig=$label)=~s/$sigPattern/$1/;if ($seenLen{$sig}<=length($seq)){print $label,@content unless $out{$sig}++;};}sub Single_Label_Reader {my $file=shift;my $label_pattern=shift;die "Error:$file unreachable.\n" unless (-s $file);open IN,q{<},$file or die "Error:$file unreadable.\n";my $temp1=<IN>;my $temp2=$temp1;my @temp;return sub {return q{} if (eof(IN));while(!eof(IN)){$temp1=$temp2;@temp=();while(<IN>){if (/$label_pattern/){$temp2=$_;last;}push @temp,$_;}return [$temp1,@temp];}}}' AllExon.fa '^>\S+\s+gene=(\S+)$' > AllExon.fa.maxLen
# 注意，修改后的命令，最长的转录本如果只包括1个序列，那么只保留第一个

################
#   表格类数据，以某一列为keys组成的Group中仅保留其对应某属性（另一列）中值最大的一类 
#################### 
# 输入的示例数据 以tab分隔
Bol013455comp23252_c01001234003721605149326002279
Bol013455comp23252_c11002430012431113533.00E-125449
Bol013456comp8286_c01003750061098437510693
Bol013456comp8286_c1100544009651508544101005
Bol013456comp200_c0100344001495183813440636
Bol013456comp200_c1100322001830215132215.00E-169595
Bol013456comp345067_c01003050017748130511.00E-159564
Bol013456comp823_c010015700222623822731172.00E-77291
# 使用的脚本
perl -e '$file=shift;$key=shift;$value=shift;open IN,$file;while(<IN>){chomp;@F=split /\t/,$_;$seen{$F[$key]}=$F[$value] if ((!$seen{$F[$key]}) || $seen{$F[$key]}<$F[$value]);}seek IN,0,0;while(<IN>){chomp;@F=split /\t/,$_;print qq{$_\n} if ($F[$value]>=$seen{$F[$key]})}'  <tab_sep_file> key_cols value_cols
# 示例命令
perl -e '$file=shift;$key=shift;$value=shift;open IN,$file;while(<IN>){chomp;@F=split /\t/,$_;$seen{$F[$key]}=$F[$value] if ((!$seen{$F[$key]}) || $seen{$F[$key]}<$F[$value]);}seek IN,0,0;while(<IN>){chomp;@F=split /\t/,$_;print qq{$_\n} if ($F[$value]>=$seen{$F[$key]})}' <tab_sep_file> 0 -1
# 输出
Bol013455comp23252_c01001234003721605149326002279
Bol013456comp8286_c1100544009651508544101

# ######
#   小文件行随机化
#######
cat infile|perl -e 'BEGIN{use List::Util qw/shuffle/;}print shuffle <>' > outfile

###################
# 打印匹配行及其前'指定数目'行 相当于 linux下的 grep -B 1 -P '\bCDS\b'
##################
perl -ne 'BEGIN{$BackwordN=1;$pattern=qr/\bCDS\b/}chomp;push @a,$_;if (@a>$BackwordN+1){shift @a};print +join ("\n",(@a)),"\n" if ($_=~$pattern)' s2_gff.gb
######################
# 打印匹配行及其后'指定数目'行 grep -A 1 -P '\bCDS\b'
#########################
perl -ne 'BEGIN{$ForseeN=1;$pattern=qr/\bCDS\b/}if ($remain>0){print;$remain--};if ($_=~$pattern){print;$remain=$ForseeN}' s2_gff.gb


############
#    多个文件区别对待
###############
# -n似乎应该是while(<ARGV>) 而不是 while(<STDIN>) 这两者还是有点区别的
# 以下可以逐行遍历，或许有用，比如写一个很长的命令，要一次处理3个文件之间的映射关系....
perl -ne '$file++ if (eof(ARGV));}print $file;{'

#########################
#   按照列名提取文件多列  BUG 已修复
###########################
# 修改@colNames即可
perl -F'\t' -lane 'BEGIN{@colNames=(qq{geneID},qq{FPKM})}if($.==1){for $index(0..$#F){push @need,$index if grep {$F[$index]=~$_} @colNames}}print join qq{\t},@F[@need]' annotation.xls

####################3
#   批量提取多个序列多个区段
####################
# 批量提取多个记录多个区段的序列，如果起始坐标大于终止坐标，默认取反向互补
# 输入 表格文件 最终区段序列ID\t目标序列ID\t起始坐标\t终止坐标
Y67    Contig\1    14039    12901
9S    Contig\1    9981    9482
76    Contig\2    10743    10528
A99    Contig\3    5511    5003
# 输入fa文件
perl -e '$info=shift;$fa=shift;open INFO,$info;while(<INFO>){chomp;($printID,$id,$start,$end)=split /\s+/,$_;if($start>$end){$dire=qq{-};($start,$end)=($end,$start);}else{$dire=qq{+}}$id=~s/^(\S+).*/$1/;push @{$cutID{$id}},[$printID,$id,$start,$end,$dire]};;;;my $label_pattern=qq{>};my $file=$fa;my $ref=&Single_Label_Reader($file,$label_pattern);my $content;while($content=&$ref){my ($label,@content)=@$content;;;;chomp $label;$seq=join q{},@content;$seq=~s/\r?\n//g;$label=~s/^>(\S+).*/$1/;if($cutID{$label}){for(@{$cutID{$label}}){($printID,$id,$start,$end,$dire)=@{$_};
if($dire eq qq{+}){$mod=1}else{$mod=-1}$cutSeq=substr($seq,$start+$mod,$end-$start+1);warn qq{$start,$end-$start\t}.length($cutSeq).qq{\n};if($dire eq qq{-}){$cutSeq=&rev_com($cutSeq)};;;print qq{>$printID\t$label\t$start\t$end\t$dire\n$cutSeq\n}};;;;}}sub Single_Label_Reader {my $file=shift;my $label_pattern=shift;die qq{Error:$file unreachable.\n} unless (-s $file);open IN,q{<},$file or die qq{Error:$file unreadable.\n};my $temp1=<IN>;my $temp2=$temp1;my @temp;return sub {return q{} if (eof(IN));while(!eof(IN)){$temp1=$temp2;@temp=();while(<IN>){if (/$label_pattern/){$temp2=$_;last;}push @temp,$_;}return [$temp1,@temp];}}}sub rev_com {my $seq=shift;chomp $seq;my @seq=split qq{},$seq;@seq=reverse @seq;$seq=join q{},@seq;$seq=~tr/ATCGatcg/TAGCtagc/;return $seq;}' extract.Info fish.fna > extracted.fa

# 如果要默认截取上下游个多2000bp
perl -e '$info=shift;$fa=shift;open INFO,$info;while(<INFO>){chomp;($printID,$id,$start,$end)=split /\s+/,$_;if($start>$end){$dire=qq{-};($start,$end)=($end,$start);}else{$dire=qq{+}}$id=~s/^(\S+).*/$1/;push @{$cutID{$id}},[$printID,$id,$start,$end,$dire]};;;;my $label_pattern=qq{>};my $file=$fa;my $ref=&Single_Label_Reader($file,$label_pattern);my $content;while($content=&$ref){my ($label,@content)=@$content;;;;chomp $label;$seq=join q{},@content;$seq=~s/\r?\n//g;$label=~s/^>(\S+).*/$1/;if($cutID{$label}){for(@{$cutID{$label}}){($printID,$id,$start,$end,$dire)=@{$_};if($dire eq qq{+}){$mod=1}else{$mod=-1}$cutSeq=substr($seq,$start+$mod-2000,$end-$start+1+4000);warn qq{$start,$end-$start\t}.length($cutSeq).qq{\n};if($dire eq qq{-}){$cutSeq=&rev_com($cutSeq)};;;print qq{>$printID\t$label\t$start\t$end\t$dire\n$cutSeq\n}};;;;}}sub Single_Label_Reader {my $file=shift;my $label_pattern=shift;die qq{Error:$file unreachable.\n} unless (-s $file);open IN,q{<},$file or die qq{Error:$file unreadable.\n};my $temp1=<IN>;my $temp2=$temp1;my @temp;return sub {return q{} if (eof(IN));while(!eof(IN)){$temp1=$temp2;@temp=();while(<IN>){if (/$label_pattern/){$temp2=$_;last;}push @temp,$_;}return [$temp1,@temp];}}}sub rev_com {my $seq=shift;chomp $seq;my @seq=split qq{},$seq;@seq=reverse @seq;$seq=join q{},@seq;$seq=~tr/ATCGatcg/TAGCtagc/;return $seq;}' extract.Info fish.fna > extractedUp2000Down2000.fa

##########
#    输出fasta文件每个序列对应的长度
############
perl -0076 -ane '@F=map{s/[>\r\n]//gr}@F;$id=shift @F;print $id,qq{\t},length (join q{},@F),qq{\n} if $id' in.fa

#############
#    java jar包发布前 依赖lib中的外源jar清理
############## 
#  引用了外源的jar包之后，程序会迅速加大，但其实我们并不需要那么多依赖，删除没必要的lib中的jar是第一步
#  verbose模式 输出 (应该进行完全的操作测试，方可得到最全的列表，否则会出错) 后过滤 jar 随后输出一个列表，依据此列表，删除lib文件夹下没有在这个列表中的jar包
java -verbose -jar  C:\Users\CJ\Documents\NetBeansProjects\WebDriver\dist\WebDriver.jar | perl -ne "print if /^\[/"|perl -ne "print unless  /jre1\.8/"|perl -ne "next unless /jar\]$/;s/.*\/(.*?jar)\]/$1/;print"|perl -ne "$count{$_}++}for (sort keys %count){print}{"


############
#   依据step长度输出字符串所有后kmer子串
#############
# 子函数 返回一个匿名字符串数组

# 原始版本 -- &getKmerByLenAcStep(字符串,Kmer长度,步长)
#            -- 注意，最后一个kmer如果长度低于指定Kmer长度，不返回该Kmer
sub getKmerByLenAcStep{
    my $str = shift;
    my $len = shift;
    my $stepLen = shift;
    my @returnArr;
    for (my $start=0;$start+$len<=length($str);$start+=$stepLen){
        push @returnArr,substr($str,$start,$len);
    }
    return \@returnArr;
}

# 单行版本
sub getKmerByLenAcStep{my $str = shift;my $len = shift;my $stepLen = shift;my @returnArr;for (my $start=0;$start+$len<=length($str);$start+=$stepLen){push @returnArr,substr($str,$start,$len);}return \@returnArr;}

# 实际例子
perl -e '$str=shift;print for @{&getKmerByLenAcStep($str,5,1)};sub getKmerByLenAcStep{my $str = shift;my $len = shift;my $stepLen = shift;my @returnArr;for (my $start=0;$start+$len<=length($str);$start+=$stepLen){push @returnArr,substr($str,$start,$len);}return \@returnArr}' 123456789

###############
#   基于SAM文件统计ref的每个序列的uniq counts并输出reads的uniq mapped rate统计信息（用于表达谱差异分析
##############
# 输入SAM文件  ==== -a -strata -best 对于bowtie
# 标准输出指向counts table
# 标准错误指向最终统计信息
perl -lane 'next if /^@/;next if(0x4&$F[1]);$readAndSub{qq{$F[0]}}++}for (keys %readAndSub){$total++;next if $readAndSub{$_}>1;$uniq++;}{open IN,$ARGV;while(<IN>){@F=split (/\t/,$_);$count{$F[2]}++ if $readAndSub{qq{$F[0]}}==1}for(keys %count){print qq{$_\t$count{$_}}}}{print STDERR qq{TOTAL:$total\nUNIQ:$uniq\nUNIQ_RATE:@{[$uniq/$total]}}}{' sample1.sam 1>sample1.counts 2>sample1.stats

##############
#   汇总所有counts table并进行无表达补零操作（用于表达谱差异分析
#################
head *.counts
# ------------------------  每行两列 第一列是基因ID 第二列是Counts 中间用Tab分隔
==> J0-R-stat.txt <==
Litchi_GLEAN_10043873   113
Litchi_GLEAN_10054101   157
Litchi_GLEAN_10056191   15
Litchi_GLEAN_10000953   3403
Litchi_GLEAN_10015890   1
Litchi_GLEAN_10030049   362
Litchi_GLEAN_10005957   35
Litchi_GLEAN_10024410   128
Litchi_GLEAN_10027051   106
Litchi_GLEAN_10048601   459
#-----------------------------
# 写一个命令专门用于补零
perl -e '@AllFile=@ARGV;print qq{GeneID\t};print qq{$_\t} for @AllFile;print qq{\n};while(<ARGV>){chomp;${$ARGV}{(split qq{\t},$_)[0]}=(split qq{\t},$_)[1];$uniqID{(split qq{\t},$_)[0]}++};for $id (keys %uniqID){print qq{$id\t};for(@AllFile){if(${$_}{$id}){print ${$_}{$id}}else{print qq{0}};print qq{\t}};print "\n"}' *.counts|perl -pe 's/\t$//'> AllCount.xls

############
#   保留fastq文件指定长度的读段最优子串
###############
# 输入为 in.fq finalReadLen out.fq
perl -e 'my $infq= shift;my $len = shift;my $outfq = shift;die $usage unless (-s $infq && $len && $outfq);open IN,q{<},$infq or die qq{Can_t read $infq\n};open OUT,q{>},$outfq or die qq{Can_t write into $outfq\n};while(!eof(IN)){my $id=<IN>;my $seq=<IN>;my $plus=<IN>;my $qual=<IN>;chomp ($id,$seq,$plus,$qual);if(length($seq)>$len){my $start = &extractBestFQ($qual,$len);print OUT qq{$id\n},substr($seq,$start,$len),qq{\n},qq{$plus\n},substr($qual,$start,$len),qq{\n};}else{print OUT qq{$id\n$seq\n$plus\n$qual\n};}}sub extractBestFQ{my $qual=shift;my $len = shift;my $final_start=0;my $aQ = 0;for (my $start=0;$start+$len<=length($qual);$start++){if(&averageQ(substr($qual,$start,$len))>$aQ){$aQ = &averageQ(substr($qual,$start,$len));$final_start = $start;}}return $final_start;}sub averageQ{my $qual_str=shift;my @allQ = split q{},$qual_str;my $totalQ = 0;for (@allQ){$totalQ+=(ord($_)-33);}return $totalQ/@allQ;}' sample.fq 36 sample.36.fq

##########
#   输出fasta文件每个记录的A T G C 字数统计
###########
perl -0076 -ane '@F=map{s/[>\r\n]//gr}@F;$id=shift @F;map {$aC+=($_=~tr/Aa//);$cC+=($_=~tr/Cc//);$gC+=($_=~tr/Gg//);$tC+=($_=~tr/Tt//)} @F;print qq{$id\t$aC\t$tC\t$gC\t$cC\n} if($id);' subject.fa

#########
#    合并配对的读段文件fastq 正反读段交错 
##########
# 注意，正反剬文件中的读段在位置上必须严格对应
perl -e '$fq1=shift;$fq2=shift;$mergedfq=shift;open F1,q{<},$fq1;open F2,q{<},$fq2;while(!eof(F1)){$id1=<F1>;$seq1=<F1>;$plus1=<F1>;$qual1=<F1>;$id2=<F2>;$seq2=<F2>;$plus2=<F2>;$qual2=<F2>; unless ($id1=~/\/1$/){chomp $id1;$id1.=qq{/1\n}};unless ($id2=~/\/2$/){chomp $id2;$id2.=qq{/2\n};};print $id1,$seq1,$plus1,$qual1,$id2,$seq2,$plus2,$qual2}' p1.fq p2.fq > pineapple.fq

#########
#    统计SAM文件 CIGAR的命令
#############
perl -le 'map {($value,$name)=$_=~/(\d+)(\D)/;$count{$name}+=$value} map {/(\d+\D+)/g} $ARGV[0];print qq{$_\t$count{$_}} for keys %count' 11M1D31M1I3M1D48M2I1M1D28M
# --------------
I       3
M       122
D       3
# --------------

###########
#   fasta文件去除ID行完全重复的记录
############
perl -0076 -ane '@F=map{s/[>\r\n]//gr}@F;$id=shift @F;print qq{>$id},qq{\n},@F,qq{\n} unless $seen{$id}++' plant.pep0 > plantko.uniqId.pep0000

#############
#    合并所有文件的指定列
#############
# 如下，合并所有文件的 第2列 （$wantColminusOne+1）
perl -lane 'BEGIN{$wantColminusOne=1;}$line++;push @{$cut{$line}},$F[$wantColminusOne];$line=0 if (eof(ARGV))}for(sort {$a <=> $b} keys %cut){print join qq{\t},@{$cut{$_}}}{' *.stat.xls
# 也可以用shell命令  有一个缺点，他不能再 win 下用
source <(echo paste $(ls|xargs -I {} echo \<\(cut -f2 {}\)|paste -s)) > dampIt

###########
#    根据id文件提取第二个文件中多个id匹配行
###############
perl -e '$file1=shift;$file2=shift;open F1,$file1;open F2,$file2;while(<F1>){chomp;$seen{$_}=1}while(<F2>){($id)=(split)[0];print if $seen{$id}}' ids nr  > 22.nr.annot.xls

#############
#  根据某一列的不同值将一个文件分割为多个文件
################
perl -e '$file=shift;$col=shift;$sep=shift;$sep=qq{\t} unless defined $sep;open IN,qq{<},$file;while(<IN>){chomp;@F=split /$sep/,$_;push @{$content{$F[$col]}},$_;}for(keys %content){open OUT,qq{>},$_;print OUT qq{$_\n} for @{$content{$_}};close OUT;}' gene2module.map 1

################
#    保留高表达或者去除低表达(WGCNA) 
################
# 注意 保留高表达基因的意思是，只保留在多个样本中均有高于 给定阈值的 表达量
# 过滤规则 80%的样本中基因的fpkm值必须高于 1
#       所有样本中fpkm的 CV值(变异系数) 必须高于 0.2
perl -lane 'BEGIN{$minFpkm=1;$minMinFpkmPercent=1.0;$minCV=0.2}if($.==1){print;next}shift @F;$count=0;$sum=0;$devSum=0;for(@F){$count++ if $_>=$minFpkm;$sum+=$_}next unless $count/@F>=$minMinFpkmPercent;$mean=$sum/@F;for(@F){$devSum+=($_-$mean)**2;};$var=$devSum/(@F-1);$sd=$var**(0.5);$CV=$sd/$mean;next unless $CV>=$minCV;print' genes.fpkm_table.mod.Len200 > genes.len.200.fpkm.1.cv.0.2

# 注意 去除低表达基因的意思是，在多个样本中表达量均低下，亦即只保留至少存在少数样本表达量较高
# 过滤掉 90% 样本中 少于 10 个 counts 的
# 筛选变异系数
perl -lane 'BEGIN{$minFpkm=10;$minMinFpkmPercent=0.9;$minCV=0}if($.==1){print;next}shift @F;$count=0;$sum=0;$devSum=0;for(@F){$count++ if $_<$minFpkm;$sum+=$_}next if $count/@F>=$minMinFpkmPercent;$mean=$sum/@F;for(@F){$devSum+=($_-$mean)**2;};$var=$devSum/(@F-1);$sd=$var**(0.5);$CV=$sd/$mean;next unless $CV>=$minCV;print' AllCount.xls > filteredNoExpress

##########
#    表格类数据依据第一列，加和其他所有列，去冗余
##############
## 原始文件
a       1       2       3
b       1       2       3
a       1       2       3
b       1       2       3
c       1       2       3
d       1       2       3
a       1       2       3
## 结果文件
a       3       6       9
b       2       4       6
d       1       2       3
c       1       2       3
## 
perl -F'\t' -lane '$geneName = shift @F;push @{$num{$geneName}},[@F];} for my $geneName (keys %num){$size=@{$num{$geneName}->[0]};my @sum;for my $arrRef(@{$num{$geneName}}){for(my $i=0;$i<$size;$i++){$sum[$i]+=$arrRef->[$i];}}print qq{$geneName\t},join qq{\t},@sum;}{'  input.table

###########
#    ghostz比对到nr的表格提取query2gi.table
##############
# 提取每一个hit对应的第一个gi，类似blast2go
perl -lane '$F[2]=~s/gi\|(\d+).*?$/$1/;print $F[0],"\t",$F[2]' AllGenes2nrGhost.xls|head
# 提取所有gi ，考虑了一个hit在描述中有多余1个gi号的情况
perl -lane '@geneId = map {/gi\|(\d+)/g} $_;print $F[0],qq{\t},+join qq{\t},@geneId' AllTrans.fa2nrGhost.1e_5.xls|head

##########
#    fastqReader
###############
# 处理 Fq记录数据即可 
perl -ne 'push @fq,$_;if(@fq==4){......;@fq=()}' sample.fq
# 过滤fastq数据
perl -ne 'push @fq,$_;if(@fq==4){$_=substr($fq[1],0,10);print @fq if (tr/A//)>=7||/A{6,}/;@fq=()}' sample.fq

############
#    Linux下依据 SRA run number下载SRA数据
###########
# 
perl -e '$id=shift;$run=substr($id,0,6);system qq{wget ftp://ftp-trace.ncbi.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/$run/$id/$id.sra}' SRR492199

# 更稳健的，批量下来的命令
for file in `cat SRAlist`;do perl -e '$id=shift;$prefix=substr($id,0,3);$run=substr($id,0,6);system qq{aria2c -j 20 ftp://ftp-trace.ncbi.nih.gov/sra/sra-instant/reads/ByRun/sra/$prefix/$run/$id/$id.sra}' $file;done


############
#    快速批量统计fq.gz文件行数
###############
# 其中用到了 for 循环 ，另外安装了pigz 做并行解压缩 安装包了parallel做并行行数统计
for file in *.fq.gz;do cat <(echo -ne $file\\t) <(pigz -d -c $file|parallel --pipe wc -l|awk '{s+=$1} END {print s}') >> line.counts;done

# 如果是双端数据，那么可以检验一下正反读段函数是否正常
perl -lane 'push @tmp,$F[1];if($.%2==0){if ($tmp[0]==$tmp[1]){print qq{$F[0] is equal.}}else{print STDERR qq{$F[0] seen NOT GOOD.}};@tmp=()}' line.counts  1> sample.ok 2>sample.NOT
cat sample.NOT
 
##########
#    格式化mapman结果（mercator）
##########
# 在mercator进行mapman注释之后，注释结果中是按照层次结构，所以第一列是按照层次排列的
# 每一列都有单引号，不适合自行处理
# 此外序列id会全部变成小写，需要修改
perl -F'\t' -lane '@F=map{s/'\''//gr}@F;;;$F[2]=~s/tr(.*)_i\d+$/TR$1/;next unless $F[2]=~/^TR\d+/;;;print join qq{\t},@F[2,0,1,3]' ../mercator.results.txt|perl -lne 'print unless /^\d+/||/^\s+/' > Mapman.anno

###########
#   表达量表格做行标准化
###########
perl -lane 'print and next if $.==1;$geneId=shift @F;$sd=sd(@F);$mean=mean(@F);@F=map{($_-$mean)/$sd}@F;print join qq{\t},($geneId,@F);sub sd{return &var(@_)**(1/2)}sub var{my $sum;my $mean=&mean(@_);$sum+=(($_-$mean)**2) for @_;return $sum/(@_-1)}sub mean{&sum(@_)/@_}sub sum{my $sum;$sum+=$_ for @_;return $sum}' in.fpkm.table

##########
#   基于ID列表提取表格（考虑待提取的表格中有单ID对应多行记录） 
#############
perl -lane 'if($flag){print if $need{$F[0]}}else{$need{$F[0]}++}$flag=1 if eof(ARGV);' allDiff.ids out.Hormone-related.predict.anno 

# 如果要保留第二个文件的表头
perl -lane 'if($flag){print and $open=0 if $open;print if $seen{$F[0]}}else{$seen{$F[0]}++;$flag=$open=1 if eof ARGV;}' allDiff.ids counts.final.AverageSorted.txt

##########
#    文件批量重命名（提供一个重命名列表） 
############
# 具体重命名列表为Tab分割
# cat rename.list
A360-T01        SL
A360-T02        ML
A360-T03        YL
A360-T16        LV
A360-T17        Pe
A360-T04        YB
# 运行以下命令
perl -lane 'system qq{rename "s/$F[0]/$F[1]/" *}' rename.list

# 一个更好的命令
perl -F'\t' -lane '$rename{$F[0]}=$F[1];END{while(<*>){$file=$_;map {if($file=~/\b$_\b/){system qq{mv $file\t@{[$file=~s/$_/$rename{$_}/r]}}}} keys %rename}}' rename.map

#########
#    perl中使用shell变量
#############
export num=100
perl -le 'print $ENV{num}'


#######
#    perl批量添加fasta文件前缀 （用于多个样本分开组装后合并并用于去冗余等操作）
########
for file in T01 T02 T03 T16 T17 T04 T05 T18 T06 T07 T08;do export file;perl -pe 's/^>(\S+)\s+.*$/>$ENV{file}_$1/' ../${file}.fasta;done > Zhonghua.fa


##########
#   对表达量表格或者counts表格 依据平均值进行排序 用于查看高表达的基因和低表达的基因都是什么
###########
perl -F'\t' -lane 'print and next if $.==1;my @New=@F;push @allLine,\@New;}print join qq{\t},@{$_} for sort {@a=@{$a};@b=@{$b};shift @a;shift @b;return &mean(@a)<=>&mean(@b)} @allLine;;;sub mean{&sum(@_)/@_}sub sum{my $sum;$sum+=$_ for @_;return $sum}{' counts-0.3.txt > counts-0.3.AverageSorted.txt

###########
#    将一个文件按列分割为多个文件（第一列为ID，其余所有列列名）
##############
perl -e '$file=shift;open IN,$file;$header=<IN>;@header=split /\t/,$header;chomp @header;;$id=shift @header;$count=1;open $fileHandle[$count++],qq{>},$_ for @header;while(<IN>){chomp;@F=split /\t/,$_;$geneId=shift @F;$count=1;for(@F){$curFh = $fileHandle[$count];print $curFh qq{$geneId\t$F[$count-1]\n};$count++}}' Unigene.fpkm.txt

###########
#    双联表计算卡方值
###############
perl -e '($a,$b,$c,$d)=@ARGV;print qq{chi_square:\t},($a+$b+$c+$d)*(($a*$d-$b*$c)**2)/(($a+$b)*($c+$d)*($a+$c)*($b+$d)),qq{\n}' 29 7 9 28
------------------------
        存活    死亡
A组    29       7
B组    9        28
----------------------

########
#    整理bowtie的比对结果
##############
# bowtie比对结束会输出比对信息，杂乱，需要整理
# reads processed: 17986232
# reads with at least one reported alignment: 12614849 (70.14%)
# reads that failed to align: 5371383 (29.86%)
Reported 21954854 paired-end alignments to 1 output stream(s)

# 批量整理所有样本的比对结果
perl -lane 'BEGIN{print qq{SampleName\tTotalReads\tMappedReads\tMappedRate\tUnmappedRead\tUnmappedRate}}$num=$.%4;if($num==1){m/.*\s+(\d+)\s*/;$totalRead=$1;}if($num==2){m/.*\s(\d+)\s*\(([.0-9%]+)\)\s*$/;($mappedRead,$mappedRate)=($1,$2)}if($num==3){m/.*\s(\d+)\s*\(([.0-9%]+)\)\s*$/;($unmappedRead,$unmappedRate)=($1,$2)}if($num==0){print join qq{\t},($ARGV,$totalRead,$mappedRead,$mappedRate,$unmappedRead,$unmappedRate)}' *.bowtie.log > bowtie.mapping.stat.xls
# 上面这个命令可能有点小问题，如果遇到 -m 的话
perl -lane 'BEGIN{print qq{SampleName\tTotalReads\tMappedReads\tMappedRate\tUnmappedRead\tUnmappedRate}};$count++;$count=0 if eof(ARGV);($count) = (0) and next if $count>=4;$num=$count%4;if($num==1){m/.*\s+(\d+)\s*/;$totalRead=$1;}if($num==2){m/.*\s(\d+)\s*\(([.0-9%]+)\)\s*$/;($mappedRead,$mappedRate)=($1,$2)}if($num==3){m/.*\s(\d+)\s*\(([.0-9%]+)\)\s*$/;($unmappedRead,$unmappedRate)=($1,$2)}if($num==0){print join qq{\t},($ARGV,$totalRead,$mappedRead,$mappedRate,$unmappedRead,$unmappedRate)}' *.bowtie.log

###########
#    基于给定列名顺序调整表格列顺序 
#############
# 按照指定样本顺序对列进行排序
# 原始顺序
# GeneIdSLMLYLYBMBFBCSRoLVPeBSS
# 排序后的顺序
# GeneIdFBYLMLSLLVPeYBMBBSSCSRo
perl -F'\t' -lane 'BEGIN{@definedRank=(GeneId,FB,YL,ML,SL,LV,Pe,YB,MB,BSS,CS,Ro)}if($.==1){for(@definedRank){$rank{$_}=$count++}for(@F){push @sortedRank,$rank{$_}}}print join qq{\t},@F[@sortedRank]' Unigenes.fpkm.mod.forKmeans > Unigenes.fpkm.mod.forKmeans.sorted


#######
#   整理GeneBank文件 （分离地点）
###########
# 下载示例文件
for gi in 76365841 22506766;do wget "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nucleotide&id=$gi&rettype=gb" -O $gi;done
# 整理并输出表格
perl -0057 -ne 's/\r?\n//;print qq{$ARGV\t$1\n} if /isolation_source="(.*)"/' *

#######
#    双列文件 整理 为 0-1 交集矩阵
########
cat tmp
#

ab

cb

db

ac

xc 

# 使用命令
perl -F'\t' -lane '$values{$F[0]}->{$F[1]}=1;$colName{$F[1]}=1}@colNames=keys %colName;print join qq{\t},(qq{ID},@colNames);for(keys %values){$curId=$_;@print=();for(@colNames){if($values{$curId}->{$_}){push @print,1}else{push @print,0}}print join qq{\t},($curId,@print)}{' tmp 

# 
IDcb
c01

a11

x10

d01
#  


######
#    整理bowtie2的比对结果
#######
#
perl -lane '$L=$count%15;$count++;$count=0 if eof(ARGV);if($L==0){m/(\d+) reads; of these:/;$info{$ARGV}[$L]=$1;}elsif($L==1){m/(\d[0-9. ()%]+) were paired; of these:/;$info{$ARGV}[$L]=$1;}elsif($L==3 || $L==4 || $L==2 || $L==7 || $L==11 || $L==12 || $L==13){m/(\d[0-9. ()%]+) aligned/;$info{$ARGV}[$L]=$1}elsif($L==6 || $L==9){m/(\d[0-9. ()%]+) pairs aligned /;$info{$ARGV}[$L]=$1}elsif($L==10){m/(\d[0-9. ()%]+) mates/;$info{$ARGV}[$L]=$1}elsif($L==14){m/(\d[0-9. ()%]+) overall/;$info{$ARGV}[$L]=$1}END{print qq{SampleName\tTotalReads\tTotalReadsInArePaired\tUnConcordantMappedPair\tConcordantMappedPairOneTime\tConcordantMappedPairMoreThanOneTime\tUnConcordantMappedPaired\tDiscordantMappedPair\tUnmappedPairs\tUnmappedReadsInPairs\tUnmappedRead\tMappedReadOneTime\tMappedReadMoreThanOneTime\tTotalMappedRate};for(sort keys %info){print join qq{\t},($_,@{$info{$_}}[0,1,2,3,4,6,7,9,10,11,12,13,14])}}' T*.bowtie2.error > bowtie2.stat.xls

############
     cat -n T01.bowtie2.error
     119840419 reads; of these:

     2  19840419 (100.00%) were paired; of these:

     3    1577336 (7.95%) aligned concordantly 0 times

     4    7149099 (36.03%) aligned concordantly exactly 1 time

     5    11113984 (56.02%) aligned concordantly >1 times

     6    ----

     7    1577336 pairs aligned concordantly 0 times; of these:

     8      585397 (37.11%) aligned discordantly 1 time

     9    ----

    10    991939 pairs aligned 0 times concordantly or discordantly; of these:

    11      1983878 mates make up the pairs; of these:

    12        670139 (33.78%) aligned 0 times

    13        69663 (3.51%) aligned exactly 1 time

    14        1244076 (62.71%) aligned >1 times

    15 98.31% overall alignment rate


###########
#  整理fastqc结果，提取所有样本的读段数
###########
# 
fastqc -t 30 *.gz
for file in *.zip;do unzip "$file";done
for file in *.zip;do perl -ne 'print qq{$ARGV\t$1\n} and exit if(/Total\s*Sequences\s*(\d+)/)' ${file%.zip}/fastqc_data.txt;done > fq.readCounts.summary.txt
# 

############
#   整理STAR比对结果
##############
perl -F'\s*\|\s*' -lane 'next if /^\s*$/;$info{$ARGV}->{$F[0]}=$F[1];push @sampleRank,$ARGV unless $seen{$ARGV}++;push @colRank,$F[0] unless $Seen{$F[0]}++;END{print join qq{\t},(qq{SampleName},map {s/^\s*//r} @colRank);for $sample (@sampleRank){@print=();push @print,$sample;for (@colRank){push @print,$info{$sample}->{$_}}print join qq{\t},@print};}' T*_Log.final.out > STAR.stat.xls



# 不知道这个可以做啥
cat okok
1,1,2,3,3,4
2,4,5,5
# 
 perl -F',' -lane '$num{$_}->[$.-1]++ for @F;}for(sort {$a<=>$b} keys %num){while(@{$num{$_}}<$.){push @{$num{$_}},0}print join qq{\t},($_,map {if($_ eq q{}){0}else{$_}} @{$num{$_}})}{' okok
1       2       0
2       1       1
3       2       0
4       1       1
5       0       2


#########
#    统计fastq文件的碱基数以及Q20和Q30的碱基数
###########
perl -lne 'next unless $.%4==0;for(split //,$_){$qual=ord($_)-33;$count++;$q20++ if $qual>=20;$q30++ if $qual>=30}}print qq{Total Bases:$count\nQ20 Bases:$q20\nQ30 bases:$q30};{' in.fq

######
#    基于NCBI的GI列表或者Accession列表获取相应物种信息
##############
# 准备GI或者Accession列表文件
# cat GIorAccession.list
DY20267
XP_009140510
XP_006404894
CDY39686
XP_010525774
XP_002880536
CDY39213
XP_009137931
CDX72251
NP_565560
XP_002867341
XP_006412691
NP_194785
XP_010447699
XP_010432974
# 下载所有文件
for gi in $(cat GIorAccession.list|paste -sd' ');do wget "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nucleotide&id=$gi&rettype=gb" -O $gi;done
# 提取所有物种信息
perl -lane 's/\r?\n//;print qq{$ARGV\t$1} if /^\s*ORGANISM\s*(.*)\s*$/' *


##########
#        依据两个tab分割表格的第一列 合并两个表格
###########
perl -F'\t' -lane 'if($flag){if(exists $add{$F[0]}){print qq{$_\t$add{$F[0]}}}else{print qq{$_\t--}}}else{$add{$F[0]}=$F[1]}$flag=1 if eof(ARGV)' geneId2Gi.table.xls genes.fpkm_table > genes.fpkm_table.withGO

##########
#    TBtools GO 注释结果 转成 BinGO 输入文件
########### 
perl -F"\t" -lane "BEGIN{print qq{(species=Own)(type=All)(curator=GO)}}$_ =~ s/GO:// and print qq{$F[0] = $_} for split /,/,$F[1]" C:\Users\CJ\Desktop\浩哥转录组分析结果\GOenrichment\query2gi.table.xls > C:\Users\CJ\Desktop\浩哥转录组分析结果\GOenrichment\query2gi.table.forBinGO.xls 

##########
#    字符串去冗余（压缩相邻所有相同字符为1个，并输出个数）
#########
# 
perl -le 'print scalar split //,qq{AAAAAABBBBBCCLLLLAAAAA}=~s/(\w)\1+/$1/gr;'

########
#    等分Fasta文件为指定个数
################
# 等分Fasta文件, 指定最终文件个数，最终文件前缀
perl -0076 -ane 'BEGIN{$fileNum=10;$prefix=qq{FaSplit};for(0..$fileNum-1){open $aa[$_],q{>},qq{$prefix.$_};}}$fh=$aa[$.%$fileNum];chomp;print $fh qq{>$_} if (split /\r?\n/,$_,2)[0]' in.fa

##########
# 等分表格文件，按行等分为指定个数文件 
###############
perl -lane 'BEGIN{$fileNum=10;$prefix=qq{FaSplit};for(0..$fileNum-1){open $aa[$_],q{>},qq{$prefix.$_}}}$fh=$aa[$.%$fileNum];print $fh $_' GYG12h.collasped.aln

############
#    Fastq文件，对reads进行trimming 
##########
perl -lpe '$_=substr($_,0,100) unless $.%2' in.fq > out.fq
# 配合parallel实现并行
 
###########
#    格式话Fasta序列为 序列每行60个碱基 
############
# 最后一个是去除空格
perl -0076 -ne 'chomp;($id,@seq)=split /\r?\n/,$_,2;print qq{>$id\n},(join q{},map {s/\s+//gr} @seq)=~s/(.{60})/$1\n/gr,qq{\n}' whole.clean.fasta|grep -v -P '^\s*$'

#############
#    Fasta ID简化及信息提取
###########
# 通用性更好的简化ID的命令
perl -lpe '/(^>[^|.\s]+)/ and $_=$1'  
# header简化
perl -pe 's/^(>\S+).*?$/$1/' in.fa > out.fa
# header信息提取
perl -lne 'print s/^>(\S+)\s*(.*?)$1/$1\t$2/r if /^>/' in.fa > out.anno.txt


##########
#    整理cd-hit-est结果为 gene2trans map 用于RSEM 
#########
 perl -lne 'if(/^>Cluster (\d+)/){$gene_id=qq{Unigene_$1}}else{m/>(TR\d+\|c\d+_g\d+_i\d+)/;print qq{$gene_id\t$1}}' Litchi.Unigenes.fa.clstr 

##########
#    输出 00001 到 10000
###########
perl -le 'printf "%.5d\n",$_ for 1..10000'

##########
#    表格转置
##########
perl -F'\t' -lane 'push @row,[@F];END{$j=0;for(@F){@printRow=();$i=0;for(@row){push @printRow,$row[$i++]->[$j]}print join qq{\t},@printRow;$j++}}' in.table > out.table

########
#    fasta2table
########
# 转换fasta为fasta table
perl -0076 -ane 'chomp @F;$geneId=shift @F;print qq{$geneId\t},join (qq{},@F),qq{\n}' possible.percursor > possible.percursor.table

########
#    反向互补子函数
######
perl -le '$seq="ATCGN";print map {join q{},reverse split(q//,tr/atcgATCG/tagcTAGC/r)} $seq'

######
#    查看bam文件中 deletion 分布
#######
# 只输出Deletion起始坐标
samtools view in.bam|perl -F'\t' -lane '$line=$_;($chr,$start,$cigar)=@F[2,3,5];for(map {/(\d+\D+)/g} $cigar){($count,$flag)=m/^(\d+)(\w)$/;if($flag eq q{M}||$flag eq q{I} || $flag eq q{N}){$start+=$count;}print join qq{\t},($chr,$start,$cigar,$line) if /D$/}'
# 输出整个region
samtools view in.bam|perl -F'\t' -lane '($chr,$start,$cigar)=@F[2,3,5];for(map {/(\d+\D+)/g} $cigar){($count,$flag)=m/^(\d+)(\w)$/;if($flag eq q{M}||$flag eq q{I} || $flag eq q{N}){$start+=$count;}print join qq{\t},($chr,$start) if /D$/}'

##########
#  对矩阵各列进行求和
########
perl -lane '$sum[$_]+=$F[$_] for 0..$#F;END{print join qq{\t},@sum}' 

##########
# 对gff文件进行 group // 用于基于位置去重，或者是 给 Overlap的record 分组，方便处理 
#########
# 这里使用的 genBlast 比对到基因组，得到的gff文件
perl -lane 'next unless $F[2] eq qq{transcript};($chr,$start,$end,$score,$strand,$_)=(@F[0,3,4,5,6],$_);push @pre,[$chr,$start,$end,$score,$strand,$_];END{$count=0;for $cur (sort {$a->[0] cmp $b->[0]||$a->[1]<=>$b->[1]} @pre){if(grep {$cur->[0] eq $_->[0] and $cur->[4] eq $_->[4] and !($_->[1]>$cur->[2]||$_->[2]<$cur->[1])} @{$group{$count}}){push @{$group{$count}},$cur}else{push @{$group{++$count}},$cur}}for $gc (sort keys %group){print join qq{\t},$gc,@{$_} for @{$group{$gc}}}}' Ananas_comosus.fna_1.1c_2.3_s2_tdshift0_tddis0_tcls3.0_m2_score_i0_d16_1.gff


#############
#  将Sanger测序结果合并，.seq 文件都没有header 
############
perl -le "while(<*>){open IN,$_;@seq=<IN>;close IN;chomp @seq;print qq{>$_\n},join q{},@seq}" > Merged.fa

##########
# 解析gff文件的注释信息（MSU Rice，并转换其中的 URI编码） 
###########
# 使用了正则表达式 的 gxr 后缀
perl -MURI -lane 'next unless $F[2] eq "gene";%attr=map{split /=/,$_} split /;/,$F[8];print join qq{\t},$attr{"ID"},$attr{"Note"}=~s/(%\w{2})/URI::Escape::uri_unescape("$1")/ger' RICE_genome_MSU7_anno.gff

#######
#    整理NCBI的Genome database上下载的gff文件中mRNA到protein的ID map 
#########
# pigz -dc GCF_001876935.1_Aspof.V1_rna.fna.gz > Asparagus_officinalis.transcript.fa
# pigz -dc GCF_001876935.1_Aspof.V1_protein.faa.gz > Asparagus_officinalis.protein.fa
pigz -dc GCF_001876935.1_Aspof.V1_genomic.gff.gz >  Asparagus_officinalis.gff
perl -lane 'if($F[2] eq "exon"){$F[8]=~/Parent=([^;]+).*Genbank:([^;]+)/;$mRNA{$1}=$2}elsif($F[2] eq "CDS"){$F[8]=~/Parent=([^;]+).*Genbank:([^;]+)/;$CDS{$1}=$2}END{for (sort keys %CDS){print join qq{\t},$mRNA{$_},$CDS{$_}}}' Asparagus_officinalis.gff > Asparagus_officinalis.cds2exon.map

#######
#    统计fasta文件整体的碱基分布 ATCG.... 
#######
head Asparagus_officinalis.transcript.fa|perl -lne 'next if /^>/;$count{"\L$_"}++ for split //,$_}{print join qq{\t},$_,$count{$_} for sort keys %count'

#######
#    解析KEGG Blast2KO 的文本 - KEGG htext 格式 得到注释表格
########
perl -lne 's/<a.*?>(.*?)<\/a>/$1/g;/^A\s*<b>(.*)<\/b>/?($A=$1):/^B\s*(.*)$/?($B=$1):/^C\s*(.*?)$/?($C=$1):/^D\s*(.*?)$/?($D=$1):(/^E\s*(.*)/ and print join qq{\t},($A,$B,$C,$D,$1))' EC.result.txt > EC.result.parsed

#########
#    过滤和筛选WGCNA输出网络表格的模块信息
###########
perl -F'\t' -lane 'print and next if $.==1;print if $F[2] eq "paleturquoise" or $F[2] eq "paleturquoise_Darkred"' CytoscapeInput-nodes-all.txt > CytoscapeInput-paleturquoise-darkred.txt
perl -F'\t' -lane 'print $F[0]' CytoscapeInput-paleturquoise-darkred.txt > ids
perl -F'\t' -slane 'BEGIN{open ID,$id;$hash{$_}=1 for map {chomp;$_} <ID>}print and next if $.==1;print if $hash{$F[0]} and $hash{$F[1]} and $F[2]>=0.1' -- -id=ids CytoscapeInput-edges-all.txt > CytoscapeInput-edges-paleturquoise-darkred.txt

########
#     使用ID映射列表，替换目标文件中所有目标文本
###########
perl -slane 'BEGIN{open MAP,$map;%hash=map{chomp;split /\t/,$_} <MAP>}$line=$_;$line=~s/$_/$hash{$_}/ for keys %hash;print $line' -- -map=番茄ARFID映射表.txt originalTree.txt

##########
#     简化Fasta的ID，去除ID行空格之后的无用信息
##########
perl -i -lpe '$_=$1 if /(>\S+)/' Oryza_sativa.genome.fna.exon.fa

###########
# 批量重命名目录下的文件（由于文件名的部分匹配和替换-基于一个列表）
###############
# 

#
6418    tetrad_1
6419    tetrad_2
6420    tetrad_3
6422    tetrad_heat_1
6423    tetrad_heat_2
6424    tetrad_heat_3
##
perl -lane '$rename{$F[0]}=$F[1];END{while(<*>){chomp;system qq{mv $_ @{[$_=~s/^(\d+)[.]/$rename{$1}./r]}}}}' rename.list
#


###########
#  Gff3文件中提取ID映射关系
#############
perl -F'\t' -lane 'next unless $F[2] eq "mRNA";%attr = map {split /=/,$_} split /;/,$F[8];print join qq{\t},@attr{"ID","Name"}' Ppatens_318_v3.3.gene.gff3 > ID2Name.map


#########
#    解析cd-hits-est 的聚类结果为tab (方便后续分析)
###########
# cd-hits-est常常用于对相似度较高的序列进行聚类（基于Shared Kmers），其速度快，但是输出结果更适合人为阅读，如下
>Cluster 0
0       363nt, >Ipomoea_batatas-miRN7... *
>Cluster 1
0       360nt, >Solanum_pimpinellifolium-miR6159... *
>Cluster 2
0       359nt, >Elaeis_guineensis-miRN14... *
>Cluster 3
0       355nt, >Zea_mays-miRN21b... at 1:355:1:359/+/86.11%
1       359nt, >Zea_mays-miRN117... *
# 为了更好的解析，写了一个one-liner，我觉得以后我仍然会用到，所以记录，将聚类结果整理为tab表格，第一列为聚类ID，余下的为该类对应的所有序列ID
perl -lane 'if(/^>(.*)/){$id=$1}else{/>(.*)\.\.\. /;push @{$cluster{$id}},$1}END{for(sort {@{$cluster{$b}}<=>@{$cluster{$a}}} keys %cluster){print join qq{\t},$_,@{$cluster{$_}}}}' novel.precursor.fa.cluster.clstr > novel.precursor.fa.cluster.clstr.tab 


##########
#    整理PlantCARE的预测结果 
##############
# 提交文件到PlantCARE, 通过邮箱接收结果，结果中包括一个汇总文件
cat plantCARE_output_*.tab|perl -F'\t' -lane 'print join qq{\t},$F[0],$F[3],$F[3]+length($F[2]),$F[7] if $F[7]'|perl -F'\t' -lane 'print if $F[3]=~/responsive|response|inducibility|induction|stresses/'|perl -F'\t' -lane 'if($F[3]=~/auxin/){$F[3]=qq{Auxin}}elsif(/defense and stress responsiveness/){$F[3]=qq{Defense}}elsif(/dehydration, low-temp, salt stresses/){$F[3]=qq{Dehydration}}elsif(/gibberellin/){$F[3]=qq{Gibberellin}}elsif(/light/){$F[3]=qq{Light}}elsif(/low-temperature/){$F[3]=qq{Low-temperature}}elsif(/salicylic acid/){$F[3]=qq{Salicylic acid}}elsif(/abscisic acid/){$F[3]=qq{Abscisic acid}}elsif(/anaerobic/){$F[3]=qq{Anaerobic}}elsif(/MeJA-responsiveness/){$F[3]=qq{MeJA}}elsif(/anoxic/){$F[3]=qq{Anoxic}}elsif(/drought-inducibility/){$F[3]=qq{Drought}}elsif(/wound-responsive/){$F[3]=qq{Wound}};print join qq{\t},@F'|perl -F'\t' -lane '$F[1]-=50;$F[2]+=50;print join qq{\t},@F' > promoter.responsive.motif.position

###
#    纯粹perl批量重命名
#####
perl -lane "rename qq{$F[0].pdf},qq{$F[1].pdf}" Rename.list

##################
#    对行进行Bootstrap抽取（即指定输出行数的可重复抽取）
##################
perl -slne '$ln++}open ARGV,$ARGV;map{$_,$lnDup{$_}++} int(rand($ln))+1 for 1..$line;while(<ARGV>){chomp;$cn++;print join qq{\t},$cn,$lnDup{$cn},$_ while($lnDup{$cn}--)};{' -- -line=10 $inFile

##############
#  从eggNOGmapper注释结果中提取GO注释信息
################
# 其中 out.emapper.annotation 从eggNOGmapper网站下载
perl -F'\t' -lane 'next if /^#/;print join qq{\t},@F[0,9] unless $F[9] eq qq{-}' out.emapper.annotation | perl -lane '$F[0]=~s/\.\d+//;print join qq{\t},$F[0],$_ for split /,/,$F[1]' 


保持更新.............
=====================================================================================

 

# perl one-liner tips
1. ARGV 及 $ARGV 变量
2. 短路操作符
perl -e '$flag=0 and print qq{UnConnected...\n}'  # 当 =0 赋值的时候，不会输出
perl -e '($flag)=(0) and print qq{Connected...\n}'  # 换成列表赋值 之后，可以输出... 因为两种赋值返回数值不同


XShell + Xmanager

yangkfor

关注

2
点赞
踩
10

收藏

觉得还不错? 一键收藏
0
评论
perl 一行命令 2021-09-02

日志返回日志列表生物信息入门常用perl ‘one-liner’–by CJ chen 2014-7-18 00:05 阅读(3326)赞(76)评论(16)转载(82)分享(29)复制地址收藏夹按钮收藏更多上一篇 | 已经是最后一篇开通黄钻生物信息入门常用perl ‘one-liner’–by CJ chenconvert -delay 10 -resample 10% ls|perl -lane '/CJ\d{4}\d{4}(\d{2})(\d{2})/;print if $1>8
复制链接

扫一扫