【NGS】如何用for循环批量做mapping?

Li Shen

写在前面:

当我们做NGS数据分析时,难免会遇到这种情况:比如我们处理paired-end的RNA-seq数据时,会遇到a_R1.fa和a_R2.fa两个匹配的文件,那么这个时候如果我们用诸如TopHat之类的工具进行mapping时,会将两个文件以

tophat -G <gtf file> <ref_genome> <a_R1.fa> <a_R2.fa>

 

这种命令处理。假设我有三个paired-end的样本,我可以用

<command1> && <command2> && <command3>

 

依次执行命令,也不算麻烦。但如果有几十上百个样要处理呢?于是今天下午就想了个小方法批量处理这些文件。

 

思路:

一开始我是想将R1和R2分别放入两个文件夹中。然后用正则表达式匹配。因为同一个样本的测序文件,整个文件名只有R1和R2的区别,其他都一样。比如:CCI-GFP-1_R1_001.fa, CCI-GFP-1_R2_001.fa。那么这个时候我只要匹配R1、R2前面的字符串是否相同,如果相同,则将这两个文件名调入tophat进行处理。

 

后来发现不需要那么麻烦,用shell中的basename命令即可代替正则表达式解决问题(shell脚本写得少啊...)。

 

代码:

for file1 in ./*_R1_001.fa;do
  file2="./$(basename "$file1" _R1_001.fa)_R2_001.fa"
  if [ ! -f "$file2" ]; then
    printf 'Missing file "%s"\n' "$file2"
    continue
  fi
    tophat -p 8 -G <gtf_file> <ref_genome> -o <output_dic> $file1 $file2
done

 

总结:

其实我一开始想用正则表达式来写是因为我用python的re库成功了,然后就一直想着shell里怎么写,然后emmmm....写了半天感觉太烦了。后来百度各种查最终找到了basename这个命令……所以以后还是得多练练shell。

 

353

April 25, 2019, 4:54 p.m.

Comments:

You haven't logged in yet, Please Log in or Signup.

No Comments.