選擇1：

拿到測序下機數(shù)據(jù)，開始分析。首先創(chuàng)建目錄

mkdir {raw,clean,qc,align,mutation}
ls

我是使用filezilla將公司給的原始數(shù)據(jù)上傳到raw目錄中。一共5個樣品，其原始數(shù)據(jù)存在在以下5個目錄中，結(jié)構(gòu)如下：

屏幕截圖 2021-11-25 093004.jpg

可以看到BAJ1和BAJ2樣品目錄下包含不止一對.fq.gz文件，這是因為有時候測序數(shù)據(jù)量不夠加測產(chǎn)生的，這樣樣品的數(shù)據(jù)實在不同的lane中產(chǎn)出的，所以會出現(xiàn)多對測序結(jié)果。為方便接下來步驟的分析，現(xiàn)在首先要將多個lane 的fastq.gz合并：

#首先根據(jù)公司給的每個文件夾中的MD5.txt文件校驗一下數(shù)據(jù)的完整性
cd BAJ1/
md5sum -c MD5.txt 

屏幕截圖 2021-11-25 094650.jpg

數(shù)據(jù)很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小環(huán)境，借助find命令和xargs命令，對五個目錄下的原始.gz文件進行fastqc

conda activate wes
find  ~/www/lxw/raw/* -name *.gz| xargs fastqc -t 10 -o ./

然后將所有的fastq的qc結(jié)果移動到qc文件夾

mv *fastqc* ../qc

然后trimmomatic去除接頭和低質(zhì)量reads,如果是coonda安裝的話，直接運行下面代碼：具體參照生信軟件 | Trimmomatic (測序數(shù)據(jù)質(zhì)控) - 知乎 (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每個文件的質(zhì)控后，還需要用trimmomatic進行去除質(zhì)量較差的測序結(jié)果。我覺得這個過程比較繁瑣，查了一下，其實有一個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如果實在小環(huán)境中用conda安裝的，shell腳本在寫的時候一定要加上conda activate wes這句，不然fastp不在環(huán)境變量中，運行會報錯說找不到這個命令
接下來的比對就可以基于fastp運行結(jié)束后生成的cleandata(即XXX.fastq.fq.gz文件）