Shell三剑客 与 核酸序列
2022-09-22 本文已影响0人
生信石头
写在前面：我是想着给课题组师弟师妹们培训出的一些例题。哈哈哈，没想到也难到了自己，因此花上了一些时间思考着这些平时我们用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 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'
事实上，这个方法我也想不到，全靠百度找答案。<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" > tmp
grep -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。也不需要我们折腾着造轮子，但是知识往往都是在造轮子的过程中积累。
