#设置参考基因相关文件变量,方便后续使用
REF=$refdir/Arabidopsis_thaliana.TAIR10.dna.chromosome.4.fa
BWA_INDEX=$refdir/Arabidopsis_thaliana.TAIR10.dna.chromosome.4.fa
GFF=$refdir/Arabidopsis_thaliana.TAIR10.47.chromosome.4.gff3
GTF=$refdir/Arabidopsis_thaliana.TAIR10.47.chromosome.4.gtf
INDEX_FAI=$refdir/Arabidopsis_thaliana.TAIR10.dna.chromosome.4.fa.fai
PICARD_DICT=$refdir/Arabidopsis_thaliana.TAIR10.dna.chromosome.4.dict
(base) [niujiajia@localhost ref]$ cd $workdir
mkdir 1.fastqc
fastqc $datadir/*.gz -o $workdir/1.fastqc
(base) [niujiajia@localhost work]$ cd 1.fastqc/
(base) [niujiajia@localhost 1.fastqc]$ ll
总用量 12072
-rw-rw-r--. 1 niujiajia niujiajia 609995 10月 24 22:06 p1_r1_fastqc.html
-rw-rw-r--. 1 niujiajia niujiajia 909592 10月 24 22:06 p1_r1_fastqc.zip
-rw-rw-r--. 1 niujiajia niujiajia 617020 10月 24 22:06 p1_r2_fastqc.html
-rw-rw-r--. 1 niujiajia niujiajia 917533 10月 24 22:06 p1_r2_fastqc.zip
-rw-rw-r--. 1 niujiajia niujiajia 611095 10月 24 22:06 p2_r1_fastqc.html
-rw-rw-r--. 1 niujiajia niujiajia 910913 10月 24 22:06 p2_r1_fastqc.zip
-rw-rw-r--. 1 niujiajia niujiajia 616308 10月 24 22:07 p2_r2_fastqc.html
-rw-rw-r--. 1 niujiajia niujiajia 922345 10月 24 22:07 p2_r2_fastqc.zip
-rw-rw-r--. 1 niujiajia niujiajia 610253 10月 24 22:07 pool1_r1_fastqc.html
-rw-rw-r--. 1 niujiajia niujiajia 907660 10月 24 22:07 pool1_r1_fastqc.zip
-rw-rw-r--. 1 niujiajia niujiajia 615201 10月 24 22:07 pool1_r2_fastqc.html
-rw-rw-r--. 1 niujiajia niujiajia 917174 10月 24 22:07 pool1_r2_fastqc.zip
-rw-rw-r--. 1 niujiajia niujiajia 614258 10月 24 22:07 pool2_r1_fastqc.html
-rw-rw-r--. 1 niujiajia niujiajia 953299 10月 24 22:07 pool2_r1_fastqc.zip
-rw-rw-r--. 1 niujiajia niujiajia 622044 10月 24 22:07 pool2_r2_fastqc.html
-rw-rw-r--. 1 niujiajia niujiajia 976910 10月 24 22:07 pool2_r2_fastqc.zip
#1.安装fastp
conda install -c "bioconda/label/cf201901" fastp
#2.数据质控:对原始序列进行去接头,删除低质量的reads等等
cd $workdir #回到工作目录
mkdir 2.data_qc
cd 2.data_qc
for i in p1 p2 pool1 pool2; do
echo "fastp -t 4 -q 10 -u 50 -l 10 -i $datadir/${i}_r1.fq.gz -I $datadir/${i}_r2.fq.gz -o ${i}_1.clean.fq.gz -O ${i}_2.clean.fq.gz -h ${i}.html -j ${i}.json"
fastp fastp -t 4 -q 10 -u 50 -l 10 --length_limit 0 -i $datadir/${i}_r1.fq.gz -I $datadir/${i}_r2.fq.gz -o ${i}_1.clean.fq.gz -O ${i}_2.clean.fq.gz -h ${i}.html -j ${i}.json
done
python $scriptdir/qc_stat.py -d $workdir/2.data_qc/ -o $workdir/2.data_qc/ -p all_sample_qc
(base) [niujiajia@localhost 2.data_qc]$ ll
总用量 557080
-rw-rw-r--. 1 niujiajia niujiajia 340 10月 24 23:34 all_sample_qc.tsv
-rw-rw-r--. 1 niujiajia niujiajia 71356114 10月 24 23:31 p1_1.clean.fq.gz
-rw-rw-r--. 1 niujiajia niujiajia 73151726 10月 24 23:31 p1_2.clean.fq.gz
-rw-rw-r--. 1 niujiajia niujiajia 480130 10月 24 23:31 p1.html
-rw-rw-r--. 1 niujiajia niujiajia 128043 10月 24 23:31 p1.json
-rw-rw-r--. 1 niujiajia niujiajia 71446954 10月 24 23:31 p2_1.clean.fq.gz
-rw-rw-r--. 1 niujiajia niujiajia 74259794 10月 24 23:31 p2_2.clean.fq.gz
-rw-rw-r--. 1 niujiajia niujiajia 480340 10月 24 23:31 p2.html
-rw-rw-r--. 1 niujiajia niujiajia 128242 10月 24 23:31 p2.json
-rw-rw-r--. 1 niujiajia niujiajia 67870872 10月 24 23:31 pool1_1.clean.fq.gz
-rw-rw-r--. 1 niujiajia niujiajia 69695559 10月 24 23:31 pool1_2.clean.fq.gz
-rw-rw-r--. 1 niujiajia niujiajia 479972 10月 24 23:31 pool1.html
-rw-rw-r--. 1 niujiajia niujiajia 127838 10月 24 23:31 pool1.json
-rw-rw-r--. 1 niujiajia niujiajia 67750448 10月 24 23:31 pool2_1.clean.fq.gz
-rw-rw-r--. 1 niujiajia niujiajia 72439232 10月 24 23:31 pool2_2.clean.fq.gz
-rw-rw-r--. 1 niujiajia niujiajia 480596 10月 24 23:31 pool2.html
-rw-rw-r--. 1 niujiajia niujiajia 128548 10月 24 23:31 pool2.json
fastp功能
usage: fastp [options] ...
options:
-i, --in1 read1 input file name (string [=])
-o, --out1 read1 output file name (string [=])
-I, --in2 read2 input file name (string [=])
-O, --out2 read2 output file name (string [=])
-6, --phred64 indicate the input is using phred64 scoring (it'll be converted to phred33, so the output will still be phred33)
-z, --compression compression level for gzip output (1 ~ 9). 1 is fastest, 9 is smallest, default is 4. (int [=4])
--stdin input from STDIN. If the STDIN is interleaved paired-end FASTQ, please also add --interleaved_in.
--stdout stream passing-filters reads to STDOUT. This option will result in interleaved FASTQ output for paired-end input. Disabled by defaut.
--interleaved_in indicate that <in1> is an interleaved FASTQ which contains both read1 and read2. Disabled by defaut.
--reads_to_process specify how many reads/pairs to be processed. Default 0 means process all reads. (int [=0])
--dont_overwrite don't overwrite existing files. Overwritting is allowed by default.
-V, --verbose output verbose log information (i.e. when every 1M reads are processed).
-A, --disable_adapter_trimming adapter trimming is enabled by default. If this option is specified, adapter trimming is disabled
-a, --adapter_sequence the adapter for read1. For SE data, if not specified, the adapter will be auto-detected. For PE data, this is used if R1/R2 are found not overlapped. (string [=auto])
--adapter_sequence_r2 the adapter for read2 (PE data only). This is used if R1/R2 are found not overlapped. If not specified, it will be the same as <adapter_sequence> (string [=])
-f, --trim_front1 trimming how many bases in front for read1, default is 0 (int [=0])
-t, --trim_tail1 trimming how many bases in tail for read1, default is 0 (int [=0])
-F, --trim_front2 trimming how many bases in front for read2. If it's not specified, it will follow read1's settings (int [=0])
-T, --trim_tail2 trimming how many bases in tail for read2. If it's not specified, it will follow read1's settings (int [=0])
-g, --trim_poly_g force polyG tail trimming, by default trimming is automatically enabled for Illumina NextSeq/NovaSeq data
--poly_g_min_len the minimum length to detect polyG in the read tail. 10 by default. (int [=10])
-G, --disable_trim_poly_g disable polyG tail trimming, by default trimming is automatically enabled for Illumina NextSeq/NovaSeq data
-x, --trim_poly_x enable polyX trimming in 3' ends.
--poly_x_min_len the minimum length to detect polyX in the read tail. 10 by default. (int [=10])
-5, --cut_by_quality5 enable per read cutting by quality in front (5'), default is disabled (WARNING: this will interfere deduplication for both PE/SE data)
-3, --cut_by_quality3 enable per read cutting by quality in tail (3'), default is disabled (WARNING: this will interfere deduplication for SE data)
-W, --cut_window_size the size of the sliding window for sliding window trimming, default is 4 (int [=4])
-M, --cut_mean_quality the bases in the sliding window with mean quality below cutting_quality will be cut, default is Q20 (int [=20])
-Q, --disable_quality_filtering quality filtering is enabled by default. If this option is specified, quality filtering is disabled
-q, --qualified_quality_phred the quality value that a base is qualified. Default 15 means phred quality >=Q15 is qualified. (int [=15])
-u, --unqualified_percent_limit how many percents of bases are allowed to be unqualified (0~100). Default 40 means 40% (int [=40])
-n, --n_base_limit if one read's number of N base is >n_base_limit, then this read/pair is discarded. Default is 5 (int [=5])
-L, --disable_length_filtering length filtering is enabled by default. If this option is specified, length filtering is disabled
-l, --length_required reads shorter than length_required will be discarded, default is 15. (int [=15])
--length_limit reads longer than length_limit will be discarded, default 0 means no limitation. (int [=0])
-y, --low_complexity_filter enable low complexity filter. The complexity is defined as the percentage of base that is different from its next base (base[i] != base[i+1]).
-Y, --complexity_threshold the threshold for low complexity filter (0~100). Default is 30, which means 30% complexity is required. (int [=30])
--filter_by_index1 specify a file contains a list of barcodes of index1 to be filtered out, one barcode per line (string [=])
--filter_by_index2 specify a file contains a list of barcodes of index2 to be filtered out, one barcode per line (string [=])
--filter_by_index_threshold the allowed difference of index barcode for index filtering, default 0 means completely identical. (int [=0])
-c, --correction enable base correction in overlapped regions (only for PE data), default is disabled
--overlap_len_require the minimum length of the overlapped region for overlap analysis based adapter trimming and correction. 30 by default. (int [=30])
--overlap_diff_limit the maximum difference of the overlapped region for overlap analysis based adapter trimming and correction. 5 by default. (int [=5])
-U, --umi enable unique molecular identifer (UMI) preprocessing
--umi_loc specify the location of UMI, can be (index1/index2/read1/read2/per_index/per_read, default is none (string [=])
--umi_len if the UMI is in read1/read2, its length should be provided (int [=0])
--umi_prefix if specified, an underline will be used to connect prefix and UMI (i.e. prefix=UMI, UMI=AATTCG, final=UMI_AATTCG). No prefix by default (string [=])
--umi_skip if the UMI is in read1/read2, fastp can skip several bases following UMI, default is 0 (int [=0])
-p, --overrepresentation_analysis enable overrepresented sequence analysis.
-P, --overrepresentation_sampling one in (--overrepresentation_sampling) reads will be computed for overrepresentation analysis (1~10000), smaller is slower, default is 20. (int [=20])
-j, --json the json format report file name (string [=fastp.json])
-h, --html the html format report file name (string [=fastp.html])
-R, --report_title should be quoted with ' or ", default is "fastp report" (string [=fastp report])
-w, --thread worker thread number, default is 2 (int [=2])
-s, --split split output by limiting total split file number with this option (2~999), a sequential number prefix will be added to output name ( 0001.out.fq, 0002.out.fq...), disabled by default (int [=0])
-S, --split_by_lines split output by limiting lines of each file with this option(>=1000), a sequential number prefix will be added to output name ( 0001.out.fq, 0002.out.fq...), disabled by default (long [=0])
-d, --split_prefix_digits the digits for the sequential number padding (1~10), default is 4, so the filename will be padded as 0001.xxx, 0 to disable padding (int [=4])
-?, --help print this message