思路:
1、由SRA文件处理为fastq文件——fastq-dump --split-e
2、fastq质控
3、mapping,包括hisat2比对后用samtools排序,获得bam文件;用samtools对bam文件构建索引
4、stringtie跑gtf文件
5、准备list文件,用getTPM.py跑表达矩阵
原始文件:153个SRA文件
目 录
2、fastq-dump命令没有--split-3功能的报错
一、从SRA到fastaq
(一)shell代码
/users/songwn/user/ZY/database/NLR/TRIDC_organ_tpm/sra_to_fastaq.sh
#!/bin/bash
path="/data/songwn/YG/database/wild_emmer_RNA/Triticum_dicoccoides_tissue"
files=$(ls $path)
for f in $path/ERR*
do
#if [ ${f:0:3} = "ERR" ];then
~/user/LPZ/software/sratoolkit.2.10.5-ubuntu64/bin/fastq-dump --split-e $f
#fi
wait
done
问题代码:
#!/bin/bash
path="/data/songwn/YG/database/wild_emmer_RNA/Triticum_dicoccoides_tissue"
files=$(ls $path)
for f in ${files}
do
if [ ${f:0:3} = "ERR" ];then
~/user/LPZ/software/sratoolkit.2.10.5-ubuntu64/bin/fastq-dump --split-e $f
fi
done
1、shell中输出的路径中的文件名
for f in $path/ERR*;
do
echo $f
done
以上命令输出的文件名包含路径
2、fastq-dump命令没有--split-3功能的报错
默认版本fastq-dump v2.10.0有问题,再安装一个老版本fastq-dump v2.9.0就好了 或者替换为--split-e
3、fastq-dump输入SRA文件时要加路径
二、fastq质控
(一)准备name文件
ls|grep "ERR*"|awk -F "_" '{print $1}' > name
(二)clean.sh代码
#!/bin/bash
dos2unix /users/songwn/user/ZY/database/NLR/TRIDC_organ_tpm/clean/name
adress1="/users/songwn/user/ZY/database/NLR/TRIDC_organ_tpm/clean"
adress2="/users/songwn/user/ZY/database/NLR/TRIDC_organ_tpm"
cat ${adress1}/name |while read line
do
#mv $adress1/${line}.gz^M $adress1/${line}.gz
~/user/LPZ/software/fastp --thread 6 -i $adress2/${line}_1.fastq -o $adress1/${line}_1.fastq -I $adress2/${line}_2.fastq -O $adress1/${line}_2.fastq -j $adress1/${line}.json -h $adress1/${line}.html
done
三、mapping跑bam文件
(一)代码
/users/songwn/user/ZY/database/NLR/TRIDC_organ_tpm/mapping/mapping.sh
#!/bin/bash
ref_index=/data/songwn/SKJ/00.ref/barley_v3_index/Morex_v3_dna
adress1=/users/songwn/user/ZY/database/NLR/TRIDC_organ_tpm/clean
adress2=/users/songwn/user/ZY/database/NLR/TRIDC_organ_tpm/mapping
cat $adress1/name | while read line
do
~/anaconda3/bin/hisat2 -p 4 -q -x $ref_index -1 $adress1/${line}_1.fastq -2 $adress1/${line}_2.fastq | ~/anaconda3/bin/samtools view -bS -q 20 | ~/anaconda3/bin/samtools sort - -o $adress2/${line}_sorted.bam
~/anaconda3/bin/samtools index $adress2/${line}_sorted.bam
done
四、stringtie跑gtf
(一)代码
/users/songwn/user/ZY/database/NLR/TRIDC_organ_tpm/gtf/stringtie.sh
#!/bin/bash
adress1="/users/songwn/user/ZY/database/NLR/TRIDC_organ_tpm/mapping"
annotion="/users/songwn/user/ZY/database/NLR/all/clean/Hv_Morex.pgsb.Jul2020.gtf"
adress2="/users/songwn/user/ZY/database/NLR/TRIDC_organ_tpm/clean"
adress3="/users/songwn/user/ZY/database/NLR/TRIDC_organ_tpm/gtf"
less $adress2/name | while read line
do
~/anaconda3/bin/stringtie $adress1/${line}_sorted.bam -G $annotion -o $adress3/${line}.gtf -A $adress3/${line}.gene.tab -p 4 -l ${line} -B -e
done
(二)准备list文件
#转到gtf所在目录
cd /users/songwn/user/ZY/database/NLR/TRIDC_organ_tpm/gtf
#提取文件名处理后做name
ls|grep "gtf"|awk -F "." '{print $1}' > name
#生成路径
sed 's/ERR/\/users\/songwn\/user\/ZY\/database\/NLR\/TRIDC_organ_tpm\/gtf\/ERR/g' name > list
#name与路径名粘贴至文件
paste name list > list.txt
#将tab键转成空格
sed -i 's/\t/ /g' list.txt
#在末尾加上文件后缀
sed -i 's/$/&\.gtf/g' list.txt
(三)getTPM.py获得表达矩阵
python2.7 getTPM.py -i list.txt