Shell三剑客 与 核酸序列

写在前面：我是想着给课题组师弟师妹们培训出的一些例题。哈哈哈，没想到也难到了自己，因此花上了一些时间思考着这些平时我们用python或者其他现有的软件实现的功能，如何用三剑客解决？<br />只能说三剑客真的比我想象中还有强大吧，我对他们了解还不够，日后加深了解。

Shell命令与生物信息相关的示例用法

• CDS文件内提取对应ID的命令- grep&sed

• 检测CDS可编码能力 - grep & awk

• fastq转fasta -awk

• 序列反向互补- awk&sed

CDS文件内提取对应ID的命令

ID：Acyan01G0000500.1 <br />file : final.gene.cds.fa

# 了解文件格式less -S  final.gene.cds.fa#格式属于一行ID一行序列

grep Acyan01G0000500.1  final.gene.cds.fa -A1

• 如果要搜索多个ID：Acyan01G0000500.1 Acyan01G0000800.1 Acyan03G0000900.1 的序列？

vi idfilegrep -f idfile  final.gene.cds.fa  -A1 grep -f idfile  final.gene.cds.fa  -A1 |grep -v "^-"

了解 -f, -v的含义 "^" "$" 的含义

#grep参数学习：-A n 指定输出搜索目标的后n行-B n 指定输出搜索目标的前n行-f <file> 以文件内容为搜索目标-v 排除搜索到的行"^" 行头符号"$" 行尾符号

思考题：对于不是一行ID+一行序列的如何提取？

sed -n '/LITCHI014205.m1/, /^>/p' Lchinesis_genome.Chr.cds| sed '$d'

事实上，这个方法我也想不到，全靠百度找答案。<br />思路：我们需要打印ID匹配到的一行到下一个基因的一行。参数 -n 只输出符合要求的行，其他行都不打印；匹配具有ID和起始为">"字符之间的行内容都实现打印，最后我们的内容是会包括了下一个基因ID行，所以删除最后一行就可以实现了。

检测CDS序列的编码性

判断ATG 起始位点

grep -v ">"  final.gene.cds.fa |  grep -v "^ATG" |head　

grep -v ">"  final.gene.cds.fa |  grep -v "^ATG"　> tmpgrep -f tmp  final.gene.cds.fa -B1 |grep ">"

如何更加优雅？

grep -f tmp  final.gene.cds.fa -B1 |grep ">" > tmp &&  grep -f tmp  final.gene.cds.fa -B1 |grep ">" |head

grep -f tmp  final.gene.cds.fa -B1 |grep ">" |grep -f -  final.gene.cds.fa -B1 |grep ">" > fail.id

如果遇到非一行名字一行序列格式的文件？
File : Lchinesis_genome.Chr.cds

grep ">" Lchinesis_genome.Chr.cds -A1|head

grep ">" Lchinesis_genome.Chr.cds -A1|grep -v ">"|grep -v "^-" | grep   -v '^ATG'

grep ">" Lchinesis_genome.Chr.cds -A1|grep -v ">"|grep -v "^-" | grep   -v '^ATG' |grep  -f - Lchinesis_genome.Chr.cds -B1 |grep ">"

判断是否满足三的倍数？

awk -F "" '{if($1~/>/){print}else{print NF}}' final.gene.cds.fa
awk -F "" '{if($1~/>/){print}else{print NF}}' final.gene.cds.fa |cut -f1 -d" " |head

#整理可看格式awk -F "" '{if($1~/>/){print}else{print NF%3}}' tmp|cut -f1 -d" " |sed ':a;N; s/\n/\t/'|sed 's/>//' | sort -k2,2nr

• 参数解析：

-F 指定分隔符，上述例子制定了分割符号为无NF 是一个变量指代【域的数量】，指代以分隔符为分割的域的总数；很复杂，如果分隔符默认为"\t"，那么NF就是列数，但是由于我们分割符已经设定了为无，所以是指代一行的总数；$1 第一域， ~ 一般是用于正则表达式上是否符合的符号 // 用于正则式的包装所以第一个命令指代的是，如果行头是以">"为符号的话，我们就直接打印；如果不是，我们就统计行的总字符数（核酸长度）；
NF%3 就是核酸长度是否能整除3，如果有余数就是非3的倍数；
最后再用sed来去替换换行符为tab分隔符；关于:a;N; 这个陌生的符号，我下面有一个图，我感觉我也说的不太清楚。
最后的最后排序查看是否有余数不为0的ID

fastq-fasta转换

#行首匹配@字符awk '{if($1~/^@/){print}}' test1_R1.fq|head

#捕获测序序列awk '{if($1~/^[ATCG]/){print}}' test1_R1.fq|head

#捕获测序序列为带@和ATCG序列awk '{if($1~/^@/){print}if($1~/^[ATCG]/){print}}' test1_R1.fq |head

#整理格式awk '{if($1~/^@/){print}if($1~/^[ATCG]/){print}}' test1_R1.fq |sed 's/@/>/'|head#重定向到新文件awk '{if($1~/^@/){print}if($1~/^[ATCG]/){print}}' test1_R1.fq |sed 's/@/>/' > test_R1.fastq

相对上面的一些例子，fastq-fasta转换就变得更加简单了。

序列反向互补配对

实现反向

• file=test.seq

awk -F "" '{if($1~/>/){print}else{{for(i=NF;i>0;i--)printf $i}printf"\n"}}' test.seq|head -n2
awk -F "" '{if($1~/>/){print}else{{for(i=NF;i>0;i--)printf $i}printf"\n"}}' test_rev.seq

同样设置分割符为无，说明逐个字符查看；用到一个awk内置判断，if(){} else{}如果以">"起始，我们直接打印，否则我们就实现反向
实现反向用到了一个内置的for()循环，i=NF 定义变量i为行数；i>0 必须大于0；i-- 实现每循环减i=i-1然后逐一输出域值（碱基）; 注意这里用到的是printf不是print。其两个区别在于，printf只输出字符，不自带换行符，而print会在字符后加换行符。当完成所有else里的内容再额外补充一次换行符，进入下一个基因；

实现互补

sed '/^[ATCG]/y/ATCG/TAGC/' test_rev.seq > test_rev_comp.seq

• 解析： 针对开头异ATCG的行进行逐一替换；

sed 匹配条件/^[ATCG]/ 满足以A/T/C/G开头的行，就可以排除了ID行；然后再进行y动作的替换，这个替换在很多教程都没怎么提及；也是一个替换动作，但跟s有点不一样；它是一个逐一替换的动作，A-T,T-A,C-G,G-C 同时实现替换，不会重复替换；进一步实现了互补动作

因为即将打算讲解的内容是三剑客，所以以上都以三剑客进行整理命令；事实上同等替换的，容易理解的命令还是很多； 反向可以用rev， 替换可以用tr

ps： 真实数据是很多已经造好的车子，如：seqkit。也不需要我们折腾着造轮子，但是知识往往都是在造轮子的过程中积累。