利用GATK对RNA-seq数据做call SNP 或 INDEL分析

生信学习小达人

于 2024-10-31 19:49:25 发布

阅读量1.1k

点赞数 9

分类专栏：分析文章标签： python 开发语言

本文链接：https://blog.youkuaiyun.com/hgz2020/article/details/143415142

版权

GATK是一个很全面的工具，不仅可以做DNA-seq的数据分析，也可以用来做RNA-seq的数据分析，此为用GATK做RNA-seq的数据分析（snps 和indels）

conda activate GATK
cd GATK/rna-seq

1.数据下载

2、质控过滤(optional)

fastqc *.gz -o ./qc

由于fastqc报告结果总体合格，未进行过滤。
利用trimmomatic过滤低质量fastq片段代码参考如下

trimmomatic PE -phred33 -trimlog logfile sample_1.fastq.gz sample_2.fastq.gz out.read_1.fq.gz \
out.trim.read_1.fq.gz out.read_2.fq.gz out.trim.read_2.fq.gz \
ILLUMINACLIP:/home/shensuo/biosoft/Trimmomatic/Trimmomatic-0.36/adapters/TruSeq3-PE.fa:2:30:10 \ 
SLsampleINGWINDOW:5:20 LEADING:5 TRAILING:5 MINLEN:50

3、star比对

此处比对与重测序数据的比对存在差别，重测序数据比对软件利用的是Bwa mem

mkdir ../bam/star_log ; cd ../bam
star_index=/home/data/gma49/GATK/ref/hg38/star-index/GRCh38_gencode_v22_CTAT_lib_Apr032020.plug-n-play/ctat_genome_lib_build_dir/ref_genome.fa.star.idx

cat > SRR.list
SRR11618713
SRR11618726
SRR11618727

cat SRR.list | while read sample
do
cd star_log
fq1=../../raw/${sample}_1.fastq.gz
fq2=../../raw/${sample}_2.fastq.gz
STAR --runThreadN  4  --genomeDir $star_index  \
--twopassMode Basic --outReadsUnmapped None --chimSegmentMin 12  \
--alignIntronMax 100000 --chimSegmentReadGapM