转自：https://mp.weixin.qq.com/s/soRzCzFSpRtGKNM2UBxM4g
写在前面：我是想着给课题组师弟师妹们培训出的一些例题。哈哈哈，没想到也难到了自己，因此花上了一些时间思考着这些平时我们用python或者其他现有的软件实现的功能，如何用三剑客解决？<br />只能说三剑客真的比我想象中还有强大吧，我对他们了解还不够，日后加深了解。
Shell命令与生物信息相关的示例用法
- CDS文件内提取对应ID的命令- grep&sed
- 检测CDS可编码能力 - grep & awk
- fastq转fasta -awk
- 序列反向互补- awk&sed
CDS文件内提取对应ID的命令
ID：Acyan01G0000500.1
file : final.gene.cds.fa
grep Acyan01G0000500.1 final.gene.cds.fa -A1
- 如果要搜索多个ID：Acyan01G0000500.1 Acyan01G0000800.1 Acyan03G0000900.1 的序列？
vi idfile
grep -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'
事实上，这个方法我也想不到，全靠百度找答案。
思路：我们需要打印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" > tmp
如何更加优雅？
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
#整理可看格式
- 参数解析：
-F 指定分隔符，上述例子制定了分割符号为无
[图片上传失败...(image-9ba726-1663904121035)]
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行；
因为即将打算讲解的内容是三剑客，所以以上都以三剑客进行整理命令；事实上同等替换的，容易理解的命令还是很多； 反向可以用rev， 替换可以用tr
ps： 真实数据是很多已经造好的车子，如：seqkit。也不需要我们折腾着造轮子，但是知识往往都是在造轮子的过程中积累。
网友评论