生成連續(xù)的自然數(shù)
seq 1 20
#輸出是1喊崖,2……20,分隔符是換行符
seq -s "," 1 20
#輸出是1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20
分割符是,
seq -s " " -w 1 13
#輸出是01 02 03 04 05 06 07 08 09 10 11 12 13
這是寬格式數(shù)字板鬓,分隔符是空格
seq -w 1 20
#輸出01,02,……20 默認(rèn)分割符是換行符咨油,輸出的數(shù)字前面是有0的,用于寬格式柒爵。
查找指定大小的文件或文件夾
find . -name "*" -type f -size 0c | xargs -n 1 rm -rf
#刪除大小為0的文件夾
find . -name "*" -size 0c | xargs -n 1 rm -rf
#刪除大小為0的文件和文件夾
find . -name "*" -type f -size +10G | xargs -n 1 rm -rf
#找到大于10G的文件,并刪除役电。
find . -name "*.fpkm |xargs -i {} cp /tmp/fpkm/"
#復(fù)制以.fpkm結(jié)尾的文件到新的目錄/temp/fpkm目錄里。
sed使用shell變量
sed "s/${filename}/${filename}.new/g"
#使用shell變量的時(shí)候,需要使用雙引號(hào)棉胀。
行尾添加字符TAIL
sed 's/$/&TAIL/g'
行首添加字符HEAD
sed 's/^/HEAD/g'
行首和行尾同時(shí)添加
sed 's/^/START/;s/$/END/g'
匹配字符串行首添加其他字符串
sed -i "/gene/ s/^/g\*/"
file.txt #file.txt文件中包含gene
的行法瑟,行首添加上g*
字符串,這是直接修改源文件唁奢。
匹配奇數(shù)行或偶數(shù)行
cat 1.txt|sed -n '1~2p'
#獲取奇數(shù)行
cat 1.txt|sed -n '2~2p'
#獲取偶數(shù)行
匹配單引號(hào)(使用單引號(hào)包括\
來(lái)轉(zhuǎn)義'
)
sed 's/^>/sed -i '\''s\//g;' cds.fa
#替換的字符中有單引號(hào)
序列大小寫轉(zhuǎn)換
sed '/^[a-z]/s/[a-z]/\u&/g' novel.fa
#轉(zhuǎn)換fasta的序列從小寫變?yōu)榇髮懽帜?br>
sed '/^[A-Z]/s/[A-Z]/\l&/g' novel.fa
#轉(zhuǎn)換大寫為小寫字母
sed匹配指定字符之間的內(nèi)容
cat kegg1.txt |sed 's/\[.*\]//g' >kegg_des
#匹配[]
中括號(hào)包含的內(nèi)容霎挟,刪除中括號(hào)和其中的字符。
cat kegg1.txt |sed 's/\(.*\)//g' >kegg_des
#匹配()
小括號(hào)包含的內(nèi)容麻掸,刪除小括號(hào)和其中的字符酥夭。
cat genome.fa|sed 's/CM[0-9]\{2,6\}//g;s/, whole genome shotgun sequence//g;' >genome.new.fa
#匹配CM后面是0-9的數(shù)字中2-6位數(shù)字,后一個(gè)是替換指定字符串, whole genome shotgun sequence
cat genome.fa|sed 's/CM[0-9]\{6\}//g;' >genome.new.fa
#匹配CM后面是0-9的數(shù)字中6位數(shù)字
獲取匹配行和前后行
grep -A 1 INT LTR.fa
#獲取包含INT的行和后一行(after)
grep -B 1 INT LTR.fa
#獲取包含INT的行和前一行(before)
grep -C 4 INT LTR.fa
#獲取包含INT的行和前后各4行(context)
去除重復(fù)行
sort -u file1
和cat file1|sort|uniq
的作用都是去除重復(fù)行脊奋,注意:uniq使用前熬北,需要sort
uniq -c file1
顯示該行出現(xiàn)的次數(shù)
uniq -u file1
僅顯示出現(xiàn)一次的行列
uniq -d file1
僅顯示重復(fù)出現(xiàn)的行
刪除空白行
cat 1.txt|tr -s '\n'
要求文件格式是unix格式。
awk字符匹配
部分匹配
awk '$1 ~/Chr/ {print $0}' genome.gff
匹配第一列包含Chr字符的gff
awk '$1 !~/Chr/ {print $0}' genome.gff
匹配第一列不包含Chr的gff
awk '$1 ~/Chr/ && $4 <= 10000000 {print $0}' genome.gff
匹配各條Chr的前10M的gff
awk '$4 >=10000000 || $4 <= 100000000 {print $0}' genome.gff
匹配前10M或者大于100M的gff
精確匹配
awk '$4 ==10000000,$4 == 100000000 {print $0}' genome.gff
匹配10M到100M的gff 注意中間的逗號(hào)表示的是范圍诚隙,,
左右的就是對(duì)應(yīng)的行的內(nèi)容讶隐。提取的是這中間的行。
awk '$1 ==Chr1,$1 == Chr3 {print $0}' genome.gff
如果文件中是按照自然的Chr1,Chr2,Chr3這樣的順序排列久又,則會(huì)提取出第1列是Chr1到Chr3的所有的行巫延。
awk if else判斷
awk '{if($4<$5)print $1,$2,$3,$4,$5,$6,"+",$8,$9;else print $1,$2,$3,$5,$4,$6,"-",$8,$9}' InGenome.gff3|tr " " "\t" >genome.gff3
查看gff文件第9列有哪些字段
awk 'BEGIN{FS=OFS="\t"} $3=="gene"{split($9, a, ";"); for(i in a){split(a[i], b, "="); if(++c[b[1]]==1) print b[1]}}' GCF_000001405.33_GRCh38.p7_genomic.gff
awk使用shell里面的變量
注意:awk使用shell變量的時(shí)候,前面用-v 變量名=值
聲明變量的賦值,后面使用的時(shí)候直接使用籽孙,不需要再使用$
烈评,awk里面普通字符串需要使用"字符串"
,變量不能被雙引號(hào)包括犯建。
abbr=$inputfile
awk -v abbr="${abbr}" '{print $1"_"abbr}' genome.id
awk輸出指定列到最后一列
分割符號(hào)是tab awk -F "\t" '{for (i=2;i<=NF;i++)printf("%s\t", $i);print ""}' file
分割符號(hào)是空格 awk -F " " '{for (i=2;i<=NF;i++)printf("%s ", $i);print ""}' file
awk里可以進(jìn)行各種算數(shù)運(yùn)算awk深入學(xué)習(xí)
awk '$3=="gene"{print $1,$4,$5,$5-$4}' test.bed
#第5列減去第4列的值
awk '$3=="gene"{print $1,$4,$5,log($5)/log(10)}' test.bed
#第5列取以10為底的對(duì)數(shù)
cat counts_num.bed| awk '{s[$1"\t"$2"\t"$3"\t"] += $4}END{ for(i in s){ print i, s[i] } }'
#計(jì)算出第1-3列行內(nèi)容一致的第4列對(duì)應(yīng)的行的和讲冠。
awk '$3=="gene"{sum+= $5-$4+1}END{print sum}' Spo.gff3
#計(jì)算所有基因的總長(zhǎng)度
awk '$3=="gene"{sum+= $5-$4+1;i+=1}END{print sum/i}' Spo.gff3
#計(jì)算所有基因的平均長(zhǎng)度
cat Spo.genome.gff3|awk '$3=="exon" {split($9,a,"="); print a[3]}'|sort|uniq -c|awk '{sum+=$1}END{print sum/NR}'
#統(tǒng)計(jì)基因的外顯子數(shù)量
cat Sp.genome.gff3|awk '$3=="gene"{split($9,a,"=");print a[2]}'
#第三列是gene,分割第9列适瓦,打印出分割后的第二個(gè)字符竿开。
paste -d "\t" <(cat Sp.genome.gff3|awk '$3=="gene"{split($9,a,"=");print a[2]}') <(cat Sp.genome.gff3|awk '$3=="gene"{print $1"\t"$4"\t"$5}') >sp.genes
#提取一個(gè)基因組的所有基因和位置信息
cat counts_num.bed| awk '{s[$1"\t"$2"\t"$3"\t"] += $4;b[$1"\t"$2"\t"$3"\t"] += 1}END{ for(i in s){ print i, s[i],b[i] } }'|awk '{print $1"\t"$2"\t"$3"\t"log($4/$5)/log(10)}'
#提取第1-3列相同的行在第4列的值的和/對(duì)應(yīng)的行數(shù)(即均值),然后再取以10為底的對(duì)數(shù)玻熙。其實(shí)$5就是基因密度否彩。
輸出第一列中字符長(zhǎng)度大于15的行
awk 'length($1)>15'
行列轉(zhuǎn)置
awk 'BEGIN{c=0;} {for(i=1;i<=NF;i++) {num[c,i] = $i;} c++;} END{ for(i=1;i<=NF;i++){str=""; for(j=0;j<NR;j++){ if(j>0){str = str" "} str= str""num[j,i]}printf("%s\n", str)} }' file
awk重命名fasta文件
awk 'FNR==NR{a[">"$1]=$2;next} $1 in a{ sub(/>/,">"a[$1]" ",$1)}1' name.list protein.fa
name.list需要有兩列,第一列是原來(lái)的ID,第二列是新的ID.protein.fa是需要重命名的fasta文件
awk只修改其中一列的內(nèi)容后輸出所有的內(nèi)容
awk -F "\t" 'OFS="\t"{if ($3==".")$3=$1"_"$2}1' AC11_L28.vcf
此處是當(dāng)?shù)谌惺?code>.時(shí)嗦随,替換為第1列下劃線第2列列荔。
awk -F "\t" 'OFS="\t"{if (NF>4)$3=$1"_"$2}1' AC11_L28.vcf
當(dāng)列數(shù)大于4的時(shí)候敬尺,替換第3列的內(nèi)容為第1列下劃線第2列。最后面的是1是個(gè)省略符贴浙,等價(jià)于 /.*
按行合并文件
cat file1 file2 >file0
按列合并文件
paste file1 file2 >file0
指定分隔符砂吞,按列合并文件
paste -d "\t" file1 file2 >file0
cut獲取指定列
cut -f 2-8 file1
獲取以tab分割的第2到8列
cut -d, -f 1-3 3.csv
獲取以,
分割的csv文件的第1到3列。
sed輸出文件指定行
sed -n '3,101p' 3.csv
獲取第3到101行崎溃。
cat 3.csv|sed '1,5d'
輸出刪除1-5行后的內(nèi)容蜻直。
cat 3.csv|sed -e '3,$d'
刪除第3到最后1行
cat 3.csv|sed '2,5c No 2-5 number'
替換第2-5行的內(nèi)容為No 2-5 number
cat 3.csv|sed '1c\test'
替換第1行的內(nèi)容為test
cat 3.csv|sed '1a\test1'
在第1行下一行插入新的一行內(nèi)容是test1
cat 3.csv|sed '1,3a\test1'
在第1-3行每行之后插入新一行
cat 3.csv|sed '1i/test1'
在第1行上一行插入新的一行內(nèi)容是test1
cat 3.csv | sed -e '/motif/d'
刪除包含字符串’motif‘的行。
sed '1~2 s/aaa/bbb/' filename
修改奇數(shù)行的內(nèi)容
sed '2~2 s/aaa/bbb/' filename
修改偶數(shù)行的內(nèi)容
sed '3~5 s/aaa/bbb/' filename
修改第3行的內(nèi)容袁串,間隔5行的第3行概而,即3,8,13,18……行囱修。(等差數(shù)列)
合并兩個(gè)文件(合并之前赎瑰,先用sort排序)
join plant.tab 2.tab -t $'\t' >test.tab
指定輸出分隔符為tab
sort -k2 cdd.xls >cdd.xls.sort
根據(jù)第2列進(jìn)行排序
sort -k1 cddid.tbl >cddid.sort
根據(jù)第1列進(jìn)行排序
join -a 1 -1 2 -2 1 cdd.xls.sort cddid.sort -t $'\t' >cdd.final.xls
匹配第1個(gè)文件的第2列(-1 2
)和第2個(gè)文件的第1列(-2 1
)一致的行,同時(shí)輸出第1個(gè)文件中未匹配的行-a 1
,使用tab作為分割符$'\t'
join -o 1.1 2.3 cdd.xls.sort cddid.sort -t $'\t' >cdd
輸出匹配的第一個(gè)文件的第一列-o 1.1
和第二個(gè)文件的第3列2.3
匹配兩個(gè)文件
grep -f file1 file2
#匹配file1和file2中包含相同內(nèi)容的行,輸出file2中這些行的內(nèi)容蔚袍。
grep -v -f file1 file2
#匹配file1和file2中不同內(nèi)容的行乡范,輸出file2中這些行的內(nèi)容。
比較兩個(gè)文件的不同
comm -1 file1 file2
文件需要先sort啤咽,按列輸出不同晋辆。參數(shù):-123 -1
表示不顯示只在file1中出現(xiàn)的列,-2
不顯示只在file2中出現(xiàn)的列宇整,-3
不顯示file1和file2共有的列瓶佳。
comm -3 file1 file2 '
文件需要先sort,輸出的是文件1和2共有的列鳞青,但是是分布在兩列窄绒,需要用sed替換tab塑猖。
diff log2013.log log2014.log -y -W 50
按行輸出兩個(gè)文件的不同慌核,-y是以并列的方式顯示文件的異同之處 -W 50是指定欄寬押逼。
檢測(cè)數(shù)據(jù)格式(檢測(cè)字段數(shù)是否一致) (datamash需要單獨(dú)安裝)
datamash check <test.tab
文件壓縮和解壓縮pigz pigz支持多線程,速度比gzip,bzip2快的多胶惰。
壓縮 pigz -k -9 -p 48 sample.bam > bam.tar.gz
-k是保留源文件 -9是指定壓縮比例傻工,-11是最大壓縮比例,-0是最小壓縮比率孵滞。 -p 指定線程數(shù)量
打包壓縮 tar cvf - 目錄名 | pigz -k -9 -p 48 > bam.tar.gz
解壓縮unpigz -p 48 bam.tar.gz
多線程(大數(shù)據(jù)量的時(shí)候中捆,能大大加快速度)
wget https://ftp.gnu.org/gnu/parallel/parallel-20190122.tar.bz2
tar jxvf parallel-20190122.tar.bz2
cd parallel-20190122
./configure -prefix=~/software/
make -j 8
make install
parallel支持標(biāo)準(zhǔn)輸出運(yùn)行和命令調(diào)用運(yùn)行泄伪。
cat jobs2run
jobs2run是要運(yùn)行的命令腳本
bzip2 oldstuff.tar
oggenc music.flac
opusenc ambiance.wav
convert bigfile.tiff small.jpeg
ffmepg -i foo.avi -v:b 12000k foo.mp4
xsltproc --output build/tmp.fo style/dm.xsl src/tmp.xml
bzip2 archive.tar
多線程方式1:parallel --jobs 6 < jobs2run
多線程方式2:
bzip原始
cat bigfile.bin | bzip2 --best > compressedfile.bz2
bzip使用多線程
cat bigfile.bin | parallel --pipe --recend '' -k bzip2 --best > compressedfile.bz2
grep原始
grep pattern bigfile.txt
grep多線程
cat bigfile.txt | parallel --block 1M --pipe grep 'pattern'
1M是指每個(gè)cpu運(yùn)行的100萬(wàn)行。會(huì)把任務(wù)按照這個(gè)值匿级,分配給各個(gè)cpu.
awk原始
cat rands20M.txt | awk '{s+=$1} END {print s}'
awk多線程
cat rands20M.txt | parallel --pipe awk \'{s+=\$1} END {print s}\' | awk '{s+=$1} END {print s}'
wc原始
wc -l bigfile.txt
wc多線程
cat bigfile.txt | parallel --pipe wc -l | awk '{s+=$1} END {print s}'
sed原始
sed s^old^new^g bigfile.txt
sed多線程
cat bigfile.txt | parallel --pipe sed s^old^new^g
速度對(duì)比
grep>awk>sed,cut
sed和cut速度比較慢蟋滴。grep -f 也比較慢染厅。sed不支持多線程,sed是按行操作脓杉,多線程會(huì)出問(wèn)題糟秘。
多行fasta轉(zhuǎn)單行
cat test.fa | tr '\n' '\t' |sed 's/\t>/\n>/g'|sed 's/\t/\n/'|sed 's/\t//g'|sed '$s/$/$\n/'
先替換\n為\t.然后把\t>這個(gè)是頭部轉(zhuǎn)回\n>,注意第二次sed是只替換一次,就把每行第一個(gè)\t換成\n,然后把全局的\t刪除球散。第4個(gè)sed是在文件末尾添加上換行符(多次使用sed,速度慢)
awk '/^>/&&NR>1{print ""}{printf "%s",/^>/?$0"\n":$0}END{print ""}' test.fa
awk速度比較快
命令講解:
/^>/&&NR>1
匹配>
開頭且行號(hào)大于1的輸出換行符
/^>/?$0"\n":$0
是判斷是否>
開頭,是輸出本行和換行散庶,否蕉堰,則輸出本行。再用printf打印出來(lái)悲龟,但是這樣輸出的每個(gè)序列的末尾沒(méi)有換行符屋讶。所以才有前一個(gè)判斷只要不是第一行>
,都先輸出一個(gè)換行符。
END{print ""}
是在文檔最后輸出一個(gè)換行须教。
更多生物信息學(xué)的one liner https://github.com/stephenturner/oneliners#etc
文本行數(shù)超過(guò)3萬(wàn)行皿渗,linux許多命令會(huì)報(bào)錯(cuò)Argument list too long
。
解決方案轻腺,使用find
和xargs
參考地址乐疆。
find ./result_dir -name "*.kaks" |xargs cat |cut -f 1,3,4,5 | grep -v 'Sequence' >All_kaks.txt
find ./result_dir -name "*.axt.one-line"|while read id;do calculate_4DTV_correction.pl $id >${id%%one-line}4dtv;done
計(jì)算4DTv。
多進(jìn)程并行運(yùn)行
多線程工具parallel是系統(tǒng)自帶的贬养,ParaFly一般的linux系統(tǒng)也自帶
ParaFly的有點(diǎn)就是挤土,可以重跑。比如運(yùn)行下面的命令误算,用于批量壓縮文件仰美,當(dāng)然pigz可以代替gzip,實(shí)現(xiàn)多線程
ls *.fq | while read fq; do echo "gzip ${fq}";done > gzip_fq.sh
nohup ParaFly -c gzip_fq.sh -CPU 8 -failed_cmds gzip_fq.sh.failed &
如果上面的壓縮命令失敗了,會(huì)輸出到gzip_fq.sh.failed文件儿礼,調(diào)整后再次nohup ParaFly -c gzip_fq.sh -CPU 8 -failed_cmds gzip_fq.sh.failed &
即可自動(dòng)跳過(guò)已經(jīng)運(yùn)行成功的命令咖杂,只運(yùn)行識(shí)別的命令。這在流程的斷點(diǎn)繼續(xù)分析中非常有用蚊夫。
nohup parallel -j 8 < 2.sh &
一次運(yùn)行腳本2.sh里的8行命令,即使用8個(gè)cpu诉字。2.sh里是有n行需要并行執(zhí)行的程序。
快速刪除空行的速度對(duì)比
genome.fa文件大小是4G.
運(yùn)行命令 | 時(shí)間 |
---|---|
time grep -v '^$' genome.fa >genome.grep.fa | 13.45s |
time sed '/^$/d' genome.fa >genome.sed.fa | 19.565s |
time tr -s "n" genome.fa >genome.tr.fa | >3min |
time awk '/./ {print}' genome.fa >genome.awk.fa | 21.971s |
time sed -i '/^$/d' genome.fa | 18.693s |