选择1:
拿到测序下机数据,开始分析。首先创建目录
mkdir {raw,clean,qc,align,mutation}
ls
我是使用filezilla将公司给的原始数据上传到raw目录中。一共5个样品,其原始数据存在在以下5个目录中,结构如下:
可以看到BAJ1和BAJ2样品目录下包含不止一对.fq.gz文件,这是因为有时候测序数据量不够加测产生的,这样样品的数据实在不同的lane中产出的,所以会出现多对测序结果。为方便接下来步骤的分析,现在首先要将多个lane 的fastq.gz合并:
#首先根据公司给的每个文件夹中的MD5.txt文件校验一下数据的完整性
cd BAJ1/
md5sum -c MD5.txt
数据很OK-,开始合并fastq.gz,即将所有的_1.fq.gz合并为一个,所有的_2.fq.gz合并为并一个
cat BAJ1_FKDN210304375-1A_H52KTDSX2_L3_1.fq.gz BAJ1_FKDN210304375-1A_HFK5MDSX2_L4_1.fq.gz >BAJ1_two_lane_L3L4_1.fq.gz
cat BAJ1_FKDN210304375-1A_H52KTDSX2_L3_2.fq.gz BAJ1_FKDN210304375-1A_HFK5MDSX2_L4_2.fq.gz >BAJ1_two_lane_L3L4_2.fq.gz
rm *FKDN* MD5.txt
同样将BAJ2样品的fastq文件也合并
然后进入wes小环境,借助find命令和xargs命令,对五个目录下的原始.gz文件进行fastqc
conda activate wes
find ~/www/lxw/raw/* -name *.gz| xargs fastqc -t 10 -o ./
然后将所有的fastq的qc结果移动到qc文件夹
mv *fastqc* ../qc
然后trimmomatic去除接头和低质量reads,如果是coonda安装的话,直接运行下面代码:具体参照生信软件 | Trimmomatic (测序数据质控) - 知乎 (zhihu.com)
trimmomatic PE file_1.fq file_2.fq output_paired_1.fq output_unpaired_1.fq output_paired_2.fq output_unpaired_2.fq ILLUMINACLIP:HiSeq-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36
选择二:fastp一步到位
上述方法fastqc每个文件的质控后,还需要用trimmomatic进行去除质量较差的测序结果。我觉得这个过程比较繁琐,查了一下,其实有一个fastp软件可一步到位,直接将rawdata处理生成cleandata,而且也能生成指控报告,代码也很简单(下面以双端测序为例)
conda install fastp
fastp -i in.R1.fq -o out.R1.fq -I in.R2.fq -O out.R2.fq
对我这次5个样品的文件,我批量fastp的代码如下(写在xx.sh脚本中):
#!/usr/bin/bash
conda activate wes
for i in `ls ~/www/lxw/raw/*/*1.fq.gz`;
do
id=${i%1.fq.gz}
fastp -t 8 -i $i -o $id"1.fastp.fq.gz" -I ${id}"2.fq.gz" -O ${id}"2.fastp.fq.gz" -h ${id}".html"
done
然后 sh xx.sh就可以了
注意:fastp如果实在小环境中用conda安装的,shell脚本在写的时候一定要加上conda activate wes这句,不然fastp不在环境变量中,运行会报错说找不到这个命令
接下来的比对就可以基于fastp运行结束后生成的cleandata(即XXX.fastq.fq.gz文件)