【linux】单行命令-第2部分：Bioinformatics
单行命令，化繁为简，重剑无锋，大巧不工。接下来学习生信实用单行命令。
基础awk和sed命令
sed: stream editor流编辑器
提取文件的2，4，5列：
awk '{print $2,$4,$5}' input.txt
输出文件第5列中等于abc123的行：
awk ` $5 == "abc123" ' input.txt
输出第5列不等于abc123的行：
awk '$5 != "abc123" ' input.txt
输出第7列以字母a-f开头的行：
awk '$7 ~ /^[a-f]/' input.txt
输出第7列中不以字母a-f开头的行：
awk '$7 !~ /^[a-f]/' input.txt
计算第2列不重复的值保存在哈希arr中（一个值只保存一次）：
awk '!arr[$2]++' input.txt
输出第3列值比第5列大的行：
awk '$3>$5' input.txt
计算文件第一列的累加值，输出最后的结果：
awk '{sum+=$1} END {print sum}' input.txt
计算第2列的平均值：
awk '{x+=$2} END {print x/NR}' input.txt
用bar替换文件中所有的foo：
sed 's/foo/bar/g' input.txt
消除行首空格或制表符：
sed 's^[ \t]*//' input.txt
消除行结尾的空格和制表符：
sed 's/[ \t]*$//' input.txt
消除行首和行尾的空格和制表符：
sed 's/^[ \t]*//;s/[ \t]*$//' input.txt
删除空行：
sed '/^$/d' input.txt
删除包含“EndOfUsefulData”的行及其后所有的行：
sed -n '/EndOfUsefulData/,$!p' input.txt
生信单行sed&awk
提取file.txt 文件Chr1中1MB和2MB的片段之间的行信息，假设第1列是Chr信息，第3列是位置信息：
# bed 文件
cat file.bed | awk '$1 == "1" ' | awk '$3>=999999' | awk '$3<=1999999'
# gff文件
cat file.gff | awk '$1 == "1" '| awk '$5 >=1000000' | awk '$5<=2000000'
统计fastq文件的一些基本信息，包括reads数、唯一的reads数、唯一reads数的比例、出现频率最多的序列及频数和所占比例：（拿到rawdata来统计时很实用）
cat myfile.fq | awk '((NR-2)%4==0{read=$1;total++;count[read]++}END{
for (read in count) {if(!max||count[read]>max){max=count[read];maxRead=read};
if(count[read]==1){unique++}};print total,unique,unique*100/total,maxRead,count[maxRead],count[[maxRead]*100/total}'
bam文件转fastq文件：
samtools view file.bam | awk 'BEGIN {FS='\t'} {print '@' $1 '\n' $10 '\n+\n' $11}' > file.fq
输出blast结果中最高得分的结果：
awk '{if (!x[$1]++) {print $0; bitscore=($14-1)} else {if($14>bitscore) print $0} }' blastout.txt
将含多条序列fasta文件分割为多个fasta（每条一个文件）：
awk '/^>/{s=++d".fa"} {print > s}' multi.fa
输出fasta文件中每条序列的序列名及其长度：
cat file.fa | awk `$0 ~ ">"{print c;c=0; printf substr($0,2,100) "\t"; } $0 !~ ">" {c+=length($0);} END { print c; }'
将fastq文件转为fasta文件：
sed -n '1~4s/^@/>/p;2~4p' file.fq > file.fa
从第2行开始，每四行取值（从FASTQ文件提取序列）：
sed -n '2~4p' file.fq
输出中剔出第1行：
awk 'NR>1' input.txt
输出第20-80行：
awk 'NR>=20&&NR<=80' input.txt
计算第2、3行列的和，并追加到每行后输出：
awk '{print $0, $2+$3}' input.txt
计算fastq文件评价read的长度
awk 'NR%4==2(sumi+=length($0)}END{print sum(NR/4)}' input.fastq
转化VSF文件为BED文件：
sed -e 's/chr//' file.vcf | awk '{OFS="\t"; if (!/^#/){print 1,1,2-1,2,2,4"/"$5,"+"}}'
sort， uniq，cut
输出带行号的内容：
cat -n file.txt
去除重复行并计数：
cat file.txt | sort -u | wc -l
找到两文件都有的行（如果两文件都是无重复行，重定向执行“wd -l”计算同样行的行数）
sort file1 file2 | uniq -d
# 安全的方法
sort -u file1 > a
sort -u file2 > b
sort a b | uniq -d
# 用comm的方法
comm -12 file1 file2
文件按照第9列数字排序（g按照常规数值，k列）：
sort -gk9 file.txt
找到第2列出现最多的字符：
cut -f2 file.txt | sort |uniq -c | sort -k1nr | head
从文件中随机选取10行：
shuf file.txt | head -n 10
输出所有的3mer DNA组合：
echo {A,C,T,G}{A,C,T,G}{A,C,T,G}
将合并的paired-end fastq文件拆分为-1和-2两个文件：（这里加上/1在/2前面）
cat interleaved.fq | paste - - - - - - - - | tee > (cut -f 1-4 | tr "\t" "\n" > deinterleaved_1.fq) |cut -f 5-8 | tr "\t" "\n" >deinterleaved_2.fq
将一个fasta文件转成一系列短的scaffolds。比如：“>Scaffold12345”，然后移除他们，保存一个去掉他们的新文件：
samtools faidx genome.fa && grep -v Scaffold genome.fa.fai | cut -f1 | xargs -n1 samtools faidx genome.fa > geome.noscaffolds.fa
显示一个隐藏的控制字符：
python -c "f = open ('file.txt' , 'r' ); f.seek(0); file = f.readlines(); print file"
find, xargs, GNU parallel
xargs：“eXtended ARGuments”的缩写，主要与find，echo和cp等命令结合使用。
find：在指定目录下查找。
GNU parallel 下载地址：https://www.gnu.org/software/parallel/
搜素文件夹及其子目录中名称为.bam的文件（目录也包括在内）：
find . -name "*.bam"
删除所有的bam文件：（谨慎操作！！！）
find . -name "*.bam" | xargs rm
将所有.txt文件重命名为.bak：（如在对*.txt做操作前用于文件备份）
find . -name "*.txt" | sed "s/\.txt$//" | xargs -i echo mv {}.txt {}.bak | sh
同时运行12个FASTQC文件：
find *fq | parallel -j 12 "fastqc {} --outdir . "
将bam文件建索引（进输出命令，并不进运行程序）：
find *.bam | parallel --dry-run 'samtools index {}'
seqktk
Seqtk专为FASTA和FASTQ而生，能快速处理FASTA或FASTQ格式序列，也可读取gzip压缩过的FASTA和FASTQ文件。下载地址：https://github.com/lh3/seqtk
将fastq转为fasta格式：
seqtk seq -a in.fq.gz > out.fa
将fastq文件中的质量值低于20的序列屏蔽掉并转为fasta格式：
# 将序列屏蔽为小写
seqtk seq -aQ64 -q20 in.fq > out.fa
# 将序列屏蔽为N
seqtk seq -aQ64 -q20 -n N in.fq > out.fa
将fasta和fastq文件格式化为每行60个字符的多行序列并去除注释信息：
seqtk seq -Cl60 in.fa > out.fa
将多行的fastq文件转为4行的fastq文件：
seqtk seq -l0 in.fq > out.fq
生成fastq或fasta的反向互补序列：
seqtk seq -r in.fq > out.fq
根据name.lst（每行一个序列名）中的序列名提取序列：
seqtk subseq in.fq name.lst > out.fq
提取reg.bed文件中所含区域的序列：
seqtk subseq in.fa reg.bed > out.fa
将reg.bed文件中所含区域的序列屏蔽为小写：
seqtk seq -M reg.bed in.fa > out.fa
使用Phred算法从两端修剪低质量的碱基：
seqtk trimfq in.fq > out.fq
从每条read的左端修剪5bp，从右端修剪10bp：
seqtk trimfq -b 5 -e 10 in.fa > out.fa
将合并的paired-end fastq文件拆分为-1和-2 两个fastq文件：
seqtk seq -l0 -1 interleaved.fq > deinterleaved_1.fq
seqtk seq -l0 -2 interleaved.fq > deinterleaved_2.fq
GFF3注释文件
输出GFF3文件所有注释的序列：
cut -s -f 1,9 yourannots.gff3 | grep $ '\t' | cut -f 1 | sort | uniq
统计GFF3文件的基因数量：
grep -c $'\tgene\t' yourannots.gff3
从GFF3文件中提取所有的基因ID：
grep $'\tgene\t' yourannots.gff3 | perl -ne 'ID=([^;]+)/ and printf("%s\n", $1)'
统计GFF3文件每个基因的长度
grep $'\tgene\t' yourannots.gff3 | cut -s -f 4,5 | perl -ne '@v = split(/\t); printf("%d\n", $v[1]-$v[0]+1)'
参考链接：
