处理野生小麦组织的表达矩阵

z_shadow_y

已于 2024-04-25 15:07:11 修改

阅读量545

点赞数

文章标签： python

于 2022-04-12 12:06:17 首次发布

本文链接：https://blog.csdn.net/z_shadow_y/article/details/123999683

版权

思路：

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文件

一、从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