转录组—上游分析_如何拿到count矩阵

## 示例如下：
├── database # 数据库存放目录，包括参考基因组，注释文件，公共数据库等
├── project  # 项目分析目录
    └── Human-16-Asthma-Trans #具体项目
        ├── data # 数据存放目录
        │   ├── cleandata # 过滤后的数据
           	│	├── trim_galore # trim_galore过滤
		   	│	└── fastp	    # fastp过滤
        │   └── rawdata # 原始数据
        ├──  Mapping # 比对目录
        │   ├── Hisat2 # Hisat比对
        │   └── Subjunc # subjunc比对
        └── Expression # 定量
            ├── featureCounts # featureCounts
            └── Salmon # salmon定量
        └── Diff_analysis # 差异分析   

wget --output-document sratoolkit.tar.gz https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/current/sratoolkit.current-ubuntu64.tar.gz
#下载速度太慢，从网页下载再上传至服务器
tar -vxzf sratoolkit.current-ubuntu64.tar.gz 
export PATH=$PATH:$PWD/sratoolkit.3.1.1-ubuntu64/bin
which fastq-dump
#/home/data/t110558/sratoolkit.3.1.1-ubuntu64/bin/fastq-dump
#测试工具包是否正常运行
fastq-dump --stdout -X 2 SRR390728

## 记得激活环境
cat SRR_download_List.txt | while read id
do
echo prefetch ${id} -O ./
done > prefetch.command

less prefetch.command
nohup bash prefetch.command &

cat sra_data/SRR_download_List.txt | while read id
do
echo "fasterq-dump -e 32 --split-files -O fq_data/ --outfile ${id}.fastq sra_data/${id}/${id}.sra"    
# 这一步使用fastq-dump或fasterq-dump都可以
echo "pigz -p 16 -f fq_data/${id}_1.fastq"
echo "pigz -p 16 -f fq_data/${id}_2.fastq"
done > sra2fq.sh

less sra2fq.sh
nohup bash sra2fq.sh &

cd /home/data/t110558/project/Human-16-Asthma-Trans/data/rawdata/qc
nohup fastqc -t 20 -o ./ ../fq_data/SRR*.fastq.gz >qc.log &
#数据整合
multiqc *.zip -o ./ -n qc_fastqc  >./multiqc.log

#样本数量少就可以直接提交到后台运行
cd $HOME/project/Human-16-Asthma-Trans/data/rawdata/fq_data
ls *gz|cut -d"_" -f 1|sort -u |while read id;do
nohup trim_galore -j 20 -q 25 --phred33 --length 36 --stringency 3 --paired -o ../../cleandata   ${id}*.gz & 
done 

ls *gz|cut -d"_" -f 1|sort -u |while read id;do
echo trim_galore -q 25 --phred33 --length 36 --stringency 3 --paired -o ../2.1.clean_fq   ${id}_*.gz 
done   > trim_galore.sh 

for i in {0..1};do ( nohup bash ../submit.sh trim_galore.sh  2 $i 1>log.trim_galore.$i.txt 2>&1 & )  ;done 

# 如果是单端 ： 
ls *gz| sort -u |while read id;do
echo trim_galore -q 25 --phred33 --length 36 --stringency 3   -o ../cleanData   ${id} 
done   > trim_galore.sh 

cat $1 |while read id
do 
	if((i%$2==$3))
	then
  $id  
  fi # check the number1 number2
  i=$((i+1))
done
# for i in {0..9};do ( nohup bash ../submit.sh trim_galore.sh  10 $i 1>log.trim_galore.$i.txt 2>&1 &)  ;done 

## 参考基因组准备:注意参考基因组版本信息
# 下载，Ensembl：http://asia.ensembl.org/index.html
# http://ftp.ensembl.org/pub/release-104/fasta/homo_sapiens/dna/

# 进入到参考基因组目录
mkdir -p $HOME/database/GRCh38.111
cd $HOME/database/GRCh38.111

# 下载基因组序列axel  curl  
nohup axel -n 100 https://ftp.ensembl.org/pub/release-111/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz >dna.log &

# 下载转录组序列
nohup axel -n 100 https://ftp.ensembl.org/pub/release-111/fasta/homo_sapiens/cdna/Homo_sapiens.GRCh38.cdna.all.fa.gz >rna.log &

# 下载基因组注释文件
nohup axel -n 100 https://ftp.ensembl.org/pub/release-111/gtf/homo_sapiens/Homo_sapiens.GRCh38.111.gtf.gz >gtf.log &

nohup axel -n 100 https://ftp.ensembl.org/pub/release-111/gff3/homo_sapiens/Homo_sapiens.GRCh38.111.gff3.gz >gff.log&

# 上述文件下载完整后，再解压；否则文件不完整就解压会报错
# 再次强调，一定要在文件下载完后再进行解压！！！
nohup gunzip Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz Homo_sapiens.GRCh38.cdna.all.fa.gz >unzip.log &

cd $HOME/database/GRCh38.111
# Hisat2构建索引，构建索引时间比较长，建议提交后台运行，一般会运行20多分钟左右
vim Hisat2Index.sh

mkdir Hisat2Index
fa=Homo_sapiens.GRCh38.dna.primary_assembly.fa
fa_baseName=GRCh38.dna
hisat2-build -p 50 -f ${fa} Hisat2Index/${fa_baseName}

# 运行
nohup bash Hisat2Index.sh >Hisat2Index.sh.log &

cd /home/data/t110558/project/Human-16-Asthma-Trans/data/cleandata

indexPrefix=/home/data/t110558/database/GRCh38.111/Hisat2Index/GRCh38.dna
ls *gz|cut -d"_" -f 1|sort -u |while read id;do
nohup  hisat2 -p 20 -x  $indexPrefix -1 ${id}*_1_val_1.fq.gz   -2 ${id}*_2_val_2.fq.gz  -S ../../Mapping/hisat2/${id}.hisat.sam 1>../../Mapping/hisat2/hisat2.log 2>&1 & 
done 

# sam文件转bam
ls *.sam|while read id ;do (nohup samtools sort -O bam -@ 20 -o $(basename ${id} ".sam").bam   ${id}   1>sam2bam.log 2>&1 & );done
# 这个过程会输出大量中间文件
rm *.sam 
# 为bam文件建立索引
ls *.bam |xargs -@ 20 -i samtools index {}

indexPrefix=/home/data/server/reference/index/hisat/mm10/genome 
ls -lh $indexPrefix*   
cat config|while read id;do
 if((i%$2==$3))
        then
  hisat2 -p 3  -x  $indexPrefix -1 ${id}*_R1_val_1.fq.gz   -2 ${id}*_combined_R2_val_2.fq.gz  -S ${id}.hisat.sam 
samtools sort -O bam -@ 3   -o ${id}.hisat.bam  ${id}.hisat.sam 
if [  ! -f  ${id}.hisat.bam    ]; then
		rm  ${id}.hisat.sam  
fi
samtools index ${id}.hisat.bam 
  fi # check the number1 number2
  i=$((i+1))
done 
  
# for i in {0..9};do ( nohup bash run_hisat2.sh t 10 $i 1>log.hisat2.$i.txt 2>&1 & )  ;done 

cd $HOME/project/Human-16-Asthma-Trans/Expression/featureCounts

## 定义输入输出文件夹
gtf=/home/data/t110558/database/GRCh38.111/Homo_sapiens.GRCh38.111.gtf
inputdir=$HOME/project/Human-16-Asthma-Trans/Mapping/hisat2/

# featureCounts对bam文件进行计数
nohup featureCounts -T 20 -p --countReadPairs -t exon -g gene_id -a $gtf -o all.id.txt $inputdir/*.bam  1>featureCounts.log 2>&1 &

# 对定量结果质控
multiqc all.id.txt.summary

# 得到表达矩阵
# 处理表头
less -S all.id.txt |grep -v '#' |cut -f 1,7- |sed 's#/home/t_rna/project/Human-16-Asthma-Trans/Mapping/Hisat2//##g' |sed 's#.Hisat_aln.sorted.bam##g' >raw_counts.txt
#修改列名
sed -i '1s#/home/data/t110558/project/Human-16-Asthma-Trans/Mapping/hisat2//##g; 1s#.hisat.bam##g' raw_counts.txt
#查看表达矩阵
head -n 100 raw_counts.txt|column -t