南京又进入一年一度的桑拿天，组里的小服务器也彻底报销了。幸好在此之前已经提前将数据迁移至单位的超算平台上。流程需要重新部署，软件需要重装，环境变量需要重新设置，还要适应超算平台上更规范的任务管理流程。超算平台要求通过qsub脚本提交任务，脚本中要设置语言环境，作业路径，调用资源以及任务名称。之前在小服务器上，直接通过nohup+命令行&提交任务，确实简单粗暴了。现已基本调制妥当，也借机重温下微生物扩增子测序的pipeline。本人生态学专业，全部围绕扩增子数据，暂无接触其他测序数据。
设置任务路径，调用资源及队列
#!/bin/sh -login
#PBS -o /home/test
#PBS -j oe
#PBS -l nodes=1:ppn=10,walltime=999:00:00,mem=30gb
#PBS -N usearch
#PBS -q cpu
cd /home/test
下面是正式的任务命令行
#新建文件夹
mkdir 01_flash/ -p
mkdir 02_cut/ -p
mkdir 03_qc -p
mkdir 04_otu -p
#flash执行序列拼接，seqtk执行反向互补，cutadapt根据引物序列识别目标序列，trimmomatic进行质控
for file in $(ls 00_raw_PE/*_R1.fq.gz | sed 's/00_raw_PE\///')
do
flash 00_raw_PE/${file%_*}_R1.fq.gz 00_raw_PE/${file%_*}_R2.fq.gz -M 200 -o ${file%_*} -d 01_flash/ -O --allow-outies >> 01_flash/${file%_*}.flash.log
seqtk seq -r 01_flash/${file%_*}.extendedFrags.fastq > 01_flash/${file%_*}.extendedFrags.rc.fastq
cat 01_flash/${file%_*}.extendedFrags.fastq 01_flash/${file%_*}.extendedFrags.rc.fastq > 01_flash/${file%_*}.fastq
cutadapt -g GTGCCAGCMGCCGCGGTAA -m 200 01_flash/${file%_*}.fastq -o 02_cut/${file%_*}.reads.cut.fastq --too-short-output reads.short.fastq --untrimmed-output 02_cut/${file%_*}.reads.no.fp.fastq > 02_cut/${file%_*}.cutadapt.report.txt
java -jar ~/programs/Trimmomatic/trimmomatic.jar SE -phred33 02_cut/${file%_*}.reads.cut.fastq 03_qc/${file%_*}.trim.fastq SLIDINGWINDOW:20:20 MINLEN:200
usearch11 -fastq_filter 03_qc/${file%_*}.trim.fastq -fastaout 03_qc/${file%_*}.fa -relabel ${file%_*}. -threads 4
done
cat 03_qc/*.fa > 04_otu/reads.fa
cd 04_otu
#使用Edgar的denoise算法进行降噪（聚类）
#一般数据量较小时使用32位的免费版usearch11
#(1) for general reads.fa
#usearch11 -fastx_uniques reads.fa -fastaout uniques.fasta -sizeout -relabel Uniq
#usearch11 -unoise3 uniques.fasta -zotus zotus.fasta -tabbedout unoise3.txt -minsize 4
#usearch11 -otutab reads.fa -zotus zotus.fasta -otutabout otu_table.txt -threads 9
#当数据量过大时，则使用64位的付费版usearch7
#(2) for reads.fa larger than 3G
usearch7 -derep_fulllength reads.fa -output uniques.fasta -sizeout -relabel Uniq
usearch11 -unoise3 uniques.fasta -zotus zotus.fasta -tabbedout unoise3.txt -minsize 4
usearch7 --usearch_global reads.fa -db zotus.fasta -strand plus -id 0.97 -uc map.uc
python ~/programs/python_scripts/uc2otutab_c.py map.uc > otu_table.txt
### 物种注释
#(1) 使用Edgar的sintax算法进行注释
#usearch11 -sintax zotus.fasta -db ~/programs/SILVA_138/SILVA_v138_16S.udb -tabbedout sintax -strand both -sintax_cutoff 0.8 -threads 9
#awk '{$2=$3=""; print $0}' sintax > zotus.sintax
#(2) 或者使用RDP classifier进行注释
java -jar ~/programs/RDPTools/classifier.jar classify -c 0.8 -t ~/programs/rdp_classifier_2.13/src/data/classifier/16srrna/rRNAClassifier.properties -f filterbyconf -h zotus.hie -o zotus.classify zotus.fasta
cat zotus.classify |sed '1d' |sed 's/""//g' |sed 's/\t/\tD_0__/;s/\t/;D_1__/2;s/\t/;D_2__/2;s/\t/;D_3__/2;s/\t/;D_4__/2;s/\t/;D_5__/2;s/\t/;D_6__/2' > zotus.classify.format
# 超算上的hdf5版本需要设置以兼容biom格式文件的处理
export HDF5_USE_FILE_LOCKING='FALSE'
#合并otu丰度表和物种注释表
biom add-metadata -i otu_table.txt --observation-metadata-fp zotus.classify.format -o otu_table_rdp.biom --sc-separated taxonomy --observation-header OTUID,taxonomy
#统计各样本的序列数
biom summarize-table -i otu_table_rdp.biom -o summary.txt
#将biome文件转换为txt文件，进行后续的统计分析可视化操作
biom convert -i otu_table_rdp.biom -o otu_table_rdp.txt --to-tsv --table-type 'OTU table' --header-key taxonomy
网友评论