1.质控
2.建立索引
3.mapping
3.1
# -p -@为线程数；-x为构建好的索引；-1 -2为双端read1和read2的fq.gz文件；|为管道符，将上一步的内容传给samtools；samtools sor是对bam文件进行排序；-o 是输出文件，可以指定路径？ ../mp_wu/${id}.sort.bam -
/home/huawei/raw_data/YSQ/RNA_seq/sort_bam/
mv P_4AL_1.clean.fq.gz 4AL_1.clean.fq.gz
mv P_4AL_2.clean.fq.gz 4AL_2.clean.fq.gz
mv P_CS_1.clean.fq.gz CS_1.clean.fq.gz
mv P_CS_2.clean.fq.gz CS_2.clean.fq.gz
ls *gz|cut -d"_" -f 1 |uniq >id.txt
INDEX="/home/huawei/raw_data/YSQ/4AL_resequence/Chinese_Spring_v1.0/Chinese_Spring_v1.0"
cat id.txt|while read id; do (hisat2 -p 8 -x ${INDEX} -1 ${id}_1.clean.fq.gz -2 ${id}_2.clean.fq.gz |samtools sort -@ 8 -o ${id}.sort.bam -) 1>./${id}.sort.bam.log 2>&1;done
#查看一下生成的bam文件
samtools samtools view 4AL_sort.bam |less -S
#将bam转回sam，用htseq-count来看基因表达
samtools view -@ 10 -h *.bam > *.sam
#featureCounts报错后，就转为HTSeq，安装在linux下用
pip install HTSeq
# -f | --format default: sam 设置输入文件的格式，该参数的值可以是sam或bam。
#-r | --order default: name 设置sam或bam文件的排序方式，该参数的值可以是name或pos。前者表示按read名进行排序，后者表示按比对的参考基因组位置进行排序。若测序数据是双末端测序，当输入sam/bam文件是按pos方式排序的时候，两端reads的比对结果在sam/bam文件中一般不是紧邻的两行，程序会将reads对的第一个比对结果放入内存，直到读取到另一端read的比对结果。因此，选择pos可能会导致程序使用较多的内存，它也适合于未排序的sam/bam文件。而pos排序则表示程序认为双末端测序的reads比对结果在紧邻的两行上，也适合于单端测序的比对结果。很多其它表达量分析软件要求输入的sam/bam文件是按pos排序的，但HTSeq推荐使用name排序，且一般比对软件的默认输出结果也是按name进行排序的。
#-s | --stranded default: yes 设置是否是链特异性测序。该参数的值可以是yes,no或reverse。no表示非链特异性测序；若是单端测序，yes表示read比对到了基因的正义链上；若是双末端测序，yes表示read1比对到了基因正义链上，read2比对到基因负义链上；reverse表示双末端测序情况下与yes值相反的结果。根据说明文件的理解，一般情况下双末端链特异性测序，该参数的值应该选择reverse（本人暂时没有测试该参数）。
#-a | --a default: 10 忽略比对质量低于此值的比对结果。在0.5.4版本以前该参数默认值是0。
#-t | --type default: exon 程序会对该指定的feature（gtf/gff文件第三列）进行表达量计算，而gtf/gff文件中其它的feature都会被忽略。此处的名词是相对自由的，建议使用符合SO惯例的名称 (sequenceontology)，如gene，repeat_region，exon， CDS等。
#-i | --idattr default: gene_id 设置feature ID是由gtf/gff文件第9列那个标签决定的；若gtf/gff文件多行具有相同的feature ID，则它们来自同一个feature，程序会计算这些features的表达量之和赋给相应的feature ID。
#-m | --mode default: union 设置表达量计算模式。该参数的值可以有union, intersection-strict and intersection-nonempty。这三种模式的选择请见上面对这3种模式的示意图。从图中可知，对于原核生物，推荐使用intersection-strict模式；对于真核生物，推荐使用union模式。
#-o | --samout 输出一个sam文件，该sam文件的比对结果中多了一个XF标签，表示该read比对到了某个feature上。
#-q | --quiet 不输出程序运行的状态信息和警告信息。
#-h | --help 输出帮助信息
htseq-count -f bam -r name -s no -a 5 -t exon -i gene_id -m intersection-nonempty /home/huawei/raw_data/YSQ/RNA_seq/sort_bam/4AL.sort.bam ./iwgsc_refseqv1.0_2017Mar13.CDS.part.gtf >counts.txt
cat id.txt|while read id; do (htseq-count -f bam -r name -s no -a 10 -t CDS -i gene_id -m intersection-nonempty /home/huawei/raw_data/YSQ/RNA_seq/sort_bam/{id}_counts.txt);done
-m 的示意图
QQ20190719-171832@2x.png
samtools的排序方式有两种（常用）
默认方式，按照染色体的位置进行排序
samtools sort test.bam default
参数-n则是根据read名进行排序。
samtools sort -n test.bam sort_left
3.2 bam 统计
ls *.bam | while read id ;do (samtools flagstat $id > ${id}_flagstat.txt);done
4.Mark PCR重复
此步的目的：
• 标记/删除PCR重复的reads
• 为后续call变异位点增加可信度，去掉假阳性
软件：GATK4(MarkDuplicates)
存在问题："-Xmx1G -Djava.io.tmpdir=./" 是做什么的？加上后占用cpu很大，去掉后报错
mkdir gatk_markdup && cd gatk_markdup
cat >gate_markdup.sh
cat id.txt|while read id;
do
gatk --java-options "-Xmx1G -Djava.io.tmpdir=./" MarkDuplicates -I ../sort_bam/${id}.sort.bam -O ${id}_marked.bam -M ${id}.metrics 1>${id}_log.mark 2>&1
done
- 将所有的reads分组
#查看帮助文档
gatk --list
# AddOrReplaceReadGroups: Assigns all the reads in a file to a single new read-group.
#查看AddOrReplaceReadGroups的子命令
gatk AddOrReplaceReadGroups
#Required Arguments （需要哪些参数，输入、输出、-LB为数据类型、-PL测序平台、-PU比对平台/软件、-SM样品名称）:
#--INPUT,-I:String Input file (BAM or SAM or a GA4GH url). Required.
#--OUTPUT,-O:File Output file (BAM or SAM). Required.
#--RGLB,-LB:String Read-Group library Required.
#--RGPL,-PL:String Read-Group platform (e.g. illumina, solid) Required.
#--RGPU,-PU:String Read-Group platform unit (eg. run barcode) Required.
#--RGSM,-SM:String Read-Group sample name Required.
cat >add_bam.sh
cat id.txt|while read id;
do
gatk --java-options "-Xmx1G -Djava.io.tmpdir=./" AddOrReplaceReadGroups -I ${id}_marked.bam -O ${id}_add.bam --LB rna -PL illumina -PU hiseq -SM ${id} 1>${id}_log.mark 2>&1
done
#samtools index 子命令建索引
cat >idex.sh
vim idex.sh
cat id.txt|while read id
do
samtools index -@ 16 ${id}_add.bam
done
mkdir vcf && cd vcf
cat >rawvcf.sh
vim rawvcf.sh
ref="/home/huawei/raw_data/YSQ/4AL_resequence/Chinese_Spring_v1.0/Chinese_Spring_v1.0.fasta"
cat id.txt|while read id
do
gatk --java-options "-Xmx1G -Djava.io.tmpdir=./" HaplotypeCaller -R $ref -I ../gatk_markdup/${id}_add.bam -O ${id}_raw.vcf 1>${id}_log.HC 2>&1
done
7.1 snp filter
cat >snp_vcf.sh
ref="/home/huawei/raw_data/YSQ/4AL_resequence/Chinese_Spring_v1.0/Chinese_Spring_v1.0.fasta"
cat id.txt|while read id
do
gatk --java-options "-Xmx8G -Djava.io.tmpdir=./" SelectVariants -R $ref -V ${id}_raw.vcf --select-type-to-include SNP -O ${id}_snp.vcf 1>${id}_log.HC 2>&1
done
ref="/home/huawei/raw_data/YSQ/4AL_resequence/Chinese_Spring_v1.0/Chinese_Spring_v1.0.fasta"
cat id.txt|while read id;do (gatk VariantFiltration -R $ref -V ${id}_snp.vcf --filter-expression "QUAL < 60.0 || QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0 || DP <10" --filter-name "Filter" -O ${id}_fsnp.vcf.gz);done
ref="/home/huawei/raw_data/YSQ/4AL_resequence/Chinese_Spring_v1.0/Chinese_Spring_v1.0.fasta"
gatk VariantFiltration -R $ref -V ${id}_snp.vcf --filter-expression "QUAL < 60.0 || QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0 || DP <10" --filter-name "Filter" -O filter.snp.vcf.gz
gatk SelectVariants --exclude-filtered true -V ${id}_fsnp.vcf.gz -O ${id}_fdsnp.vcf.gz
7.2 indel filter
cat >indel_vcf.sh
ref="/home/huawei/raw_data/YSQ/4AL_resequence/Chinese_Spring_v1.0/Chinese_Spring_v1.0.fasta"
cat id.txt|while read id
do
gatk --java-options "-Xmx8G -Djava.io.tmpdir=./" SelectVariants -R $ref -V ${id}_raw.vcf --select-type-to-include INDEL -O ${id}_indel.vcf 1>${id}_log.HC 2>&1
done
ref="/home/huawei/raw_data/YSQ/4AL_resequence/Chinese_Spring_v1.0/Chinese_Spring_v1.0.fasta"
cat id.txt|while read id; do (gatk VariantFiltration -R $ref -V ${id}_indel.vcf --filter-expression "QUAL < 60.0 || QD < 2.0 || FS > 200.0 || SOR > 10.0 || Inbreeding Coeff < -0.8 || ReadPosRankSum < -20.0" --filter-name "Filter" -O filter.indel.vcf.gz);done
gatk SelectVariants --exclude-filtered true -V filter.indel.vcf.gz -O filtered.indel.vcf.gz
计数
featureCounts -T 10 -p -t exon -g gene_id -a /home/huawei/raw_data/YSQ/RNA_seq/featureCounts/iwgsc_refseqv1.0_2017Mar13.CDS.part.gtf -o ~/all_feature.txt /home/huawei/raw_data/YSQ/RNA_seq/sort_bam/*sort.bam
网友评论