【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

参考阅读：

生信菜鸟团 - 跑BWA比对测试一下酷睿I9的CPU