數(shù)據(jù)過濾

簡介：

目的去掉低質量的堿基或接頭，一般使用兩款軟件：trim_galore和trimmomatic,他們的去除是從不合格的位置往后全部切掉。
trim_galore：它去接頭序列依賴的是Cutadapt，接頭一般出現(xiàn)在3’末端。
它參數(shù)如下：

• -q 表示設置的堿基質量閾值
• --phred33表示質量體系。phred33和phred64，目前測序平臺出的數(shù)據(jù)都是phred33.
• --length表示測序片段長度的閾值，比如設置閾值50，若去除接頭和低質量的堿基后，長度低于50bp，就直接把整條去掉。
• --stringency設置數(shù)字表示：有幾個堿基和接頭序列是有重疊，默認值為1，意思是從3`的位置檢測，出現(xiàn)一個堿基與常用接頭有重疊就把這個堿基以后的序列都去掉。
• --paird表示雙末端
• -o 輸出路徑

實際操作

• 合并文件并加路徑

raw=~/RNAseq/raw
ls `pwd`*_1* >new_1
ls `pwd`*_2* >new_2
paste new_1 new_2>conf

• 遍歷文件并進行過濾

# 首先進入腳本目錄
cd $rna/script
# 然后新建一個腳本文件，并向其中添加內容
cat >filter.sh #回憶一下前面怎么用cat新建一個腳本
# 腳本的內容是：
rna=/My_PATH/RNAseq #這里注意修改成自己的
cat $rna/raw/conf | while read i
do
    fqs=($i)
    fq1=${fqs[0]}
    fq2=${fqs[1]}
    nohup trim_galore -q 25 --phred33 --length 50 -e 0.1 --stringency 3 --paired -o $rna/clean $fq1 $fq2 &
done

bash filter.sh

sh、bash都是linux中的解釋器，但有的sh是沒有數(shù)組array解釋器的，因此當你使用sh去運行整個命令時，有可能會產生報錯：
Syntax error: "(" unexpected (expecting "done")
但這并不是說代碼寫錯了，而是選擇錯了解釋器。當然如果你的sh可以成功解釋代碼，也是可以用的

過濾結果

比對

需要數(shù)據(jù)

• 過濾并質控合格的數(shù)據(jù)
• 參考基因組、注釋文件

• 生成小數(shù)據(jù)

# 先配置clean data路徑
cd $rna/clean
ls `pwd`/*_1*gz >1.clean
ls `pwd`/*_2*gz >2.clean
paste 1.clean 2.clean >clean.conf

• 新建test文件夾，進行測試

mkdir $rna/test && cd "$_"
cat >sample_fq.sh
# 輸入以下內容
rna=/MY_PATH/rnaseq 
cat $rna/clean/clean.conf | while read i;do
    fqs=($i)
    fq1=${fqs[0]}
    fq2=${fqs[1]}
    file=`basename $fq1`
    #echo $file
    surname=${file%%_*}
    #echo $fq1 $fq2 $surname
    # 隨機選了10000條reads
    seqtk sample $fq1 10000 >test.${surname}_1.fastq
    seqtk sample $fq2 10000 >test.${surname}_2.fastq
done

basename的語法是：basename[選項][參數(shù)]
其中：
選項：為有路徑信息的文件名，如/home/test/test.txt

參數(shù)：指文件擴展名

basename用法

seqtk 是一款針對fasta/fastq 文件進行處理的小程序，有很多的功能，速度很快，很方便；
安裝conda install seqtk

子程序

seqtk seq : 用途：
　　案例1：seq 序列常規(guī)轉換

將fastq轉換成fasta：

seqtk seq -a Sample_R1.fq.gz > Sample_R1.fa

將fastq序列做反向互補分析：

seqtk seq -r Sample_R1.fq.gz > Sample_Revc_R1.fq

案例2：sample 隨機抽樣

seqtk sample -s100 Sample_R1.fq.gz 10000

可直接對壓縮文件進行序列隨機提取，在提取R1和R2兩個文件的時候，需要-s值一致，才能使提取的序列id號對應。

案例3：subseq 提取序列

根據(jù)輸入的bed文件信息，將固定區(qū)域的序列提取出來：

seqtk subseq in.fa reg.bed > out.fa

根據(jù)輸入的name list，提取相應名稱序列：

seqtk subseq in.fq name.lst > out.fq

案例4：截取序列

切除reads的前5bp，以及后10bp：

seqtk trimfq -b 5 -e 10 in.fq > out.fq

hisat2比對

mapping的目的：找到reads在基因組上的位置。

• 構建索引

mkdir $rna/align/hisat2 && cd "$_"
cat >index.sh #輸入以下內容
rna=/MY_PATH/rnaseq
genome=$rna/ref/hg19.fa
hisat2-build -p 10 $genome hg19

bash index.sh

• 比對

# 小數(shù)據(jù)的路徑
ls $rna/test/*_1* >1.test
ls $rna/test/*_2* >2.test
paste 1.test 2.test >test.conf
#合并文件

cat >hisat2_aln.sh
#創(chuàng)建索引路徑
index= $rna/align/hisat2/hg19
#讀取目標文件
cat test.conf| while read i
do
    fqs=($i)
    fq1=${fqs[0]}
    fq2=${fqs[1]}
    sample=`basename $fq1 | cut -d '_' -f1`
    hisat2 -p 10 -x $index -1 $fq1 -2 $fq2  | samtools sort -O bam \
        -@ 10  -o ${sample}_hisat.sorted.bam
    samtools index -@ 10 -b ${sample}_hisat.sorted.bam 
done

hisat2
-x 指定基因組索引
-1 指定第一個fastq文件
-2 指定第二個fastq文件
samtools
sort對bam文件進行排序。
-o 輸出文件名
index：建立索引
-b 默認下輸出是 SAM 格式文件，該參數(shù)設置輸出 BAM 格式

結果

作業(yè)

（1）
過濾后

image.png

image.png

過濾前

image.png

• hisat2
  -x 指定基因組索引
  -1 指定第一個fastq文件
  -2 指定第二個fastq文件
• samtools
  sort對bam文件進行排序。
  -p 線程數(shù)
  -o 輸出文件名
  index：建立索引
  -b 默認下輸出是 SAM 格式文件，該參數(shù)設置輸出 BAM 格式

# 小數(shù)據(jù)的路徑
ls $rna/test/*_1* >1.test
ls $rna/test/*_2* >2.test
paste 1.test 2.test >test.conf
#合并文件

cat >hisat2_aln.sh
#創(chuàng)建索引路徑
index= $rna/align/hisat2/hg19
#讀取目標文件
cat test.conf| while read i
do
    fqs=($i)
    fq1=${fqs[0]}
    fq2=${fqs[1]}
#去除文件路徑，只保留第一個文件名
    sample=`basename $fq1 | cut -d '_' -f1`
#hisat2  10線程 將1，2兩個文件比對到$index建立好索引的基因組上，運用samtools排序，輸出bam文件
    hisat2 -p 10 -x $index -1 $fq1 -2 $fq2  | samtools sort -O bam \
        -@ 10  -o ${sample}_hisat.sorted.bam
#對生成bam文件排序
    samtools index -@ 10 -b ${sample}_hisat.sorted.bam 
done

（3）

#建立文件夾
mkdir $rna/clean/trimmomatic && cd "$_" 
#運行腳本
cat >filter-2.sh
rna=/MY_PATH/rnaseq 
cat $rna/raw/conf| while read i
do
    fqs=($i)
    fq1=${fqs[0]}
    fq2=${fqs[1]}
#basename去掉目標文件絕對路徑
    tri1=`basename $fq1`
    tri2=`basename $fq2`

    nohup trimmomatic PE -phred33 \
    -trimlog trim.logfile\
    $fq1 $fq2 \
    clean.$tri1 unpaired.$tri1 \
    clean.$tri2 unpaired.$tri2 \
    SLIDINGWINDOW:5:20 LEADING:5 TRAILING:5 MINLEN:20 &
done

bash filter-2.sh

image.png