序列处理命令行笔记
2018-03-18 本文已影响26人
lxmic
1.缩短fasta格式蛋白质id行的长度
perl -lne 'if(/>/){$_=~m/>Sol.*\.[0-9]\.[0-9]/;print "$&";}else{$_=~m/[A-Z]*\*?$/;print "$&"}' mate.fa > out.fa
2.缩短fasta格式基因id行的长度
perl -lne 'if(/>/){$_=~m/>Sol.*\.[0-9]\.[0-9]/;print "$&";}else{$_=~m/[ATCG]*/;print "$&"}' mate.fa > out.fa
3.多行序列变成单行序列
perl -pe '/^>/ ? print "\n" : chomp' in.fasta | tail -n +2 > out.fasta
4.截取文件行和列
# 截取-f显示指定字段的内容，后面跟列数（1,2,3），重定向截取的内容到out文件
cut -f 1,2,3 in.file > out.file
# 截取7-20行的内容重定向到out.file
sed -n ‘7,20p’ in.file > out.file
5.排序命令sort
# 升序排列in.file
sort in.file
# 降序排列加参数-r，默认是升序
sort -r in.file
# 按数值排序，sort默认会将10<2，需要添加参数-n，才会认为10>2
sort -n in.file
# 指定列进行排序，及分隔符指定参数-k和-t
sort -k -4 -t: in.file
# 指定第四列以降序数值分隔符为tab制表符来对文件diffinfo进行排序：需要注意的是，制表符和平时的‘\t’不一样，而是$'\t'的形式
sort -nr -k 4 -t$'\t' diffinfo -o diffinfo
6.根据ID提取fasta数据集中的序列
这种方法都可以快速的拿到id匹配的序列，但是还是有点区别，主要就是perl脚本不知道是按什么顺序排列的，samtools是按idlist来逐个获取的，这一点还是要说明一下。
6.1 perl脚本
#! /usr/bin/perl -w
use strict;
die "perl $0 <lst><fa>\n" unless @ARGV==2;
my ($lst,$fa)=@ARGV;
open IN,$lst||die;
my %ha;
map{chomp;$ha{(split)[0]}=1}<IN>;
close IN;
$fa=~/gz$/?(open IN,"gzip -cd $fa|"||die):(open IN,$fa||die);
$/=">";<IN>;$/="\n";
my %out;
while(<IN>){
my $info=$1 if(/^(\S+)/);
$/=">";
my $seq=<IN>;
$/="\n";
$seq=~s/>|\r|\*//g;
print ">$info\n$seq" if(exists $ha{$info} && ! exists $out{$info});
$out{$info}=1;
}
close IN;
6.2 Samtools工具提取序列
$ samtools faidx your.fa
$ xargs samtools faidx your.fa < idfile > outfile
7.grep 命令
# 统计匹配到string的行数，选项-c
grep -c "string" file
# 只显示匹配到的内容,选项-o
grep -o “string” file
# 使用正则来匹配内容，正则需要用‘’单引号括起来
grep -o ‘[0-9]’ file
8.sed行处理命令
#截取文本特定的行,并重定向
sed -n '7,21p' file > out
#替换文本特定字符串
sed ‘s/[a-z]/[A-Z]/g’
9.批量下载番茄12条染色体测序文件
因为每个染色体都是独立的一个文件夹，而且没有一直的数值，难以用单个变量来进行循环，我还是选择了分开两部分
#下载前9条染色体序列
for ((i=1;i<=9;i++));do wget ftp://ftp.ncbi.nlm.nih.gov/genomes/Solanum_lycopersicum/CHR_0$i/les_ref_SL2.50_chr$i.fa.gz ;done
# 下载后3条染色体序列
for ((i=10;i<=12;i++)); do wget ftp://ftp.ncbi.nlm.nih.gov/genomes/Solanum_lycopersicum/CHR_$i/les_ref_SL2.50_chr$i.fa.gz; done
10. 从转录组差异表达基因中获取某个家族的geneID
# 将xlsx文件先另存为制表符间隔的文本文档，然后再处理
grep "KEYWORD" diffexpgene.txt | cut -f 1,3,4,5,7,12 | grep -v "BGI_novel" | cut -f 1 > out.file
11. 根据geneID从gff文件中获取基因的位置信息
# 根据geneid批量获取基因的染色体位置信息
grep -f geneid.file ITAG2.4_gene_models.gff3 | cut -f 1,4,5,9 | grep "ID=mRNA" | sed 's/;Name.*//g' | sed 's/ID.*://g' > gene.bed
12. 利用bedtools中的getfasta选项，根据染色体位置信息文件下载gene序列
# 下载基因序列，用于基因结构分析
bedtools getfasta -fi genome.fa -bed gene.bed -name
