【Linux 笔记】Linux 基本操作 - 03. shell脚本编程
笔记接上篇【Linux 笔记】Linux 基本操作 - 02. shell编程基础。笔记大部分源于生信技能树的B站视频教程【生信技能树】生信人应该这样学linux(更新至第14集),如有需要,可去欣赏原汁原味的视频讲解。
9. shell脚本编程
需要学会用脚本并行处理文件,解决实操中遇到的问题。处理测序数据时,可构建自己的配置文件。配置文件的一般格式为样本名 左端测序 右端测序
,可利用以下命令进行查看:
less -SN config.mapping
注意:文件命名需规范,文件命名最好不要以数字开头。
for i in `ls *.1`; do mv $i ${i%.*}; done #批量修改文件名
less -SN config.mapping #SN: 整行显示+行号
ls -lh alignment/*.bam | grep K #一般测序文件都是~G,若文件太小,则说明有问题
#样本名 左端测序 右端测序
#先取出配置文件第二列和第三列的路径对应的样本名,查看是否匹配
cut -f 1 config.mapping > 1
cut -f 2 config.mapping | while read id; do echo $(basename $id);done | cut -d"_" -f 1 >2
#basename #即取路径的最后的文件名
#对得到的文件名,用下划线分割,分割后下划线前的部分为第一列,下划线后面的部分为第二列
cut -f 3 config.mapping | while read id; do echo $(basename $id);done | cut -d"_" -f 1 >3
paste 1 2 3 | less #paste按列组合数据
grep 'NT-2' config.mapping | less -s #文件命名"-"和"_"问题;样本名字出现2次,可能是匹配错了,出现一次的没办法和别人搞混
grep 'NT-2' config.mapping | cut -f 1> wrong.sampleID #错的少的话,可以打开vim手动删除
grep 'NT-2' config.mapping | cut -f 1 | cut -d"-" -f 1 | sort | uniq -c
grep 'NT-2' config.mapping | cut -f 1 | cut -d"-" -f 1 | sort | uniq -c | grep -w 1 #将grep 'NT-2'的结果里面有运行正确的样本,需要将这些正确的样本找出来
grep 'NT-2' config.mapping | cut -f 1 | cut -d"-" -f 1 | sort | uniq -c | grep -w 1 | awk '{print $2}'>once.ID
#用找到的这些正确样本与所有样本做比较,剔除这些正确的样本,剩下的就是需要重新运行的样本了
grep -v -f once.ID wrong.sampleID | wc #用代码删除匹配错误的样本;-v反向匹配,即保留未匹配到的部分
ls *gz | grep val_1.fq.gz >fq1
ls *gz | grep val_2.fq.gz >fq2
wc fq1 fq2
paste fq1 fq2 |less #不同文件内的字符串依然按照字典序排序,因此匹配会出现错误;思考,如果在合并之前分别加入行号作为标签,然后再合并会不会出错。
#直接根据样本名构造配置文件:路径+文件名+后缀
grep 'NT-2' config.mapping | cut -f 1 > wrong.sampleID
grep -v -f once.ID wrong.sampleID > temp
mv temp wrong.sampleID
cat wrong.sampleID
cat wrong.sampleID | while read id; do ($id"\t/home/hucy/clean/$id_R1\t/home/hucy/clean/$id_R2");done
#配置文件命名变异性过大,不适合构造,还是采取分别调整fq1和fq2排序的方法
head -4 fq1 fq2
#思路:排序的内容里边不包含"_",查看sort的帮助文档,"-t"选项可以对内容进行分割后再进行排序,而"-k"选型则可以指定依据哪一列进行排序;如-k1,1是对第一列到第一列进行排序,默认按字符进行排序。
sort -tk fq2 | head
sort -t"_" -k1,1 fq1>fq1.bak
sort -t"_" -k1,1 fq2>fq2.bak
paste fq1.bak fq2.bak | less -SN
ls *_1_fq.gz|wc
ls *_1_fq.gz|sort -t"_" -k3,3 | less -S #指对第3次出现的"_"进行分割排序
ls *_1_fq.gz|sort -t"_" -k3,3 >fq1
ls *_2_fq.gz|sort -t"_" -k3,3 >fq2
paste fq1 fq2 | less -S
grep -w -f wrong.sampleID temp | grep NT-2 >c1
grep -v "NT-" temp>c2
cat c1 c2 |less -S
cat temp | while read id;do (echo $(basename $id));done| cut -d"_" -f 1>id
paste id temp >config.again #新生成的配置文件
# 重新提交命令,替换有问题的bam文件
for i in {0..7};do sbash ../step2-wes-bwa.sh ./../config.again 8 &1;done
配置文件示例
NPC10F-N /data//NPC10F-N_1.fastq.gz /data//NPC10F-N_2.fastq.gz
NPC10F-T /data//NPC10F-T_1.fastq.gz /data//NPC10F-T_2.fastq.gz
NPC15F-N /data//NPC15F-N_1.fastq.gz /data//NPC15F-N_2.fastq.gz
NPC15F-T /data//NPC15F-T_1.fastq.gz /data//NPC15F-T_2.fastq.gz
NPC29F-N /data//NPC29F-N_1.fastq.gz /data//NPC29F-N_2.fastq.gz
NPC29F-T /data//NPC29F-T_1.fastq.gz /data//NPC29F-T_2.fastq.gz
NPC34F-N /data//NPC34F-N_1.fastq.gz /data//NPC34F-N_2.fastq.gz
NPC34F-T /data//NPC34F-T_1.fastq.gz /data//NPC34F-T_2.fastq.gz
NPC37F-N /data//NPC37F-N_1.fastq.gz /data//NPC37F-N_2.fastq.gz
NPC37F-T /data//NPC37F-T_1.fastq.gz /data//NPC37F-T_2.fastq.gz
bwa比对脚本
#!/bin/bash
# usage: for i in {0..9};do (nohup bash step1-bwa.sh ./ config.fq 10 $i & );done
analysis_dir=$1
config_file=$2
number1=$3
number2=$4
INDEX=/data/bigbiosoft/GATK/resources/bundle/hg38/bwa_index/gatk_hg38
cat $config_file |while read id
do
arr=($id)
fq1=${arr[1]}
fq2=${arr[2]}
sample=${arr[0]}
if((i%$number1==$number2))
then
#####################################################
################ Step 1 : Alignment #################
#####################################################
if [ ! -f ok.step1-Alignment.$sample.status ]; then
start=$(date +%s.%N)
echo bwa $sample `date` >> $sample.log
bwa mem -t 5 -M -R "@RG\tID:$sample\tSM:$sample\tLB:WES\tPL:Illumina" $INDEX $fq1 $fq2 1>>$sample.sam 2>>/dev/null
if [ $? -eq 0 ]; then
echo "bwa succeed" $sample >> $sample.log
touch ok.step1-Alignment.$sample.status
else
echo "failed" $sample >> $sample.log
fi
echo bwa $sample `date` >> $sample.log
dur=$(echo "$(date +%s.%N) - $start" | bc)
printf "Execution time for BWA : %.6f seconds" $dur >> $sample.log
echo >> $sample.log
fi
fi
i=$((i+1))
done
脚本中涉及的小知识点,超级多,有的我也不太熟悉,自行探索吧。
批量提交用法
for i in {0..9};do (nohup bash step1-bwa.sh ./ config.fq 10 $i & );done
参考阅读: