转录组分析--hisat2+featurecounts

1 ###文件下载

单个文件
wget -c -t 0 -O SRR13107018.sra https://sra-pub-run-odp.s3.amazonaws.com/sra/SRR13107018/SRR13107018
#-c -t 配合使用可以防止下载数据的过程中链接中断的问题，-O则可以指定下载路径和文件名。

#多个文件一起下
/home/software/sratoolkit.2.9.6-1-centos_linux64/bin/prefetch -O ./ --option-file SRR_Acc_List.txt

2 ## 循环把下载的所有sra文件都变为fastq(双端测序数据)

ls *sra | while read id ;do (nohup /home/software/sratoolkit.2.9.6-1-centos_linux64/bin/fastq-dump --gzip --split-3 $id &) ;done

3 ### fastp质量控制

ls *fastq.gz | cut -d"_" -f 1 |sort -u | while read id ;do (nohup fastp -i ${id}_1.fastq.gz -I ${id}_2.fastq.gz -g -q 5 -n 15 -l 150 -u 50 -o ${id}_1.filter.fastq.gz -O ${id}_2.filter.fastq.gz -h ${id}.fastp.html &) ;done

# 去除带接头(adapter)的双端读长(默认开启，可用-A关闭)；
# 去掉单端读长中含 N(N 表示无法确定碱基信息)的碱基比例大于10%部分（默认开启）
# -q 设置碱基质量阈值，默认是15，phred质量评分≥Q15
# -n 一条read中N碱基的数量超过了多少个，那么这条read就被删除，默认是5（即这条read里有5个N）
# -g 进行PolyG尾的去除
# -u 允许百分之多少的碱基不合格（0-100），默认是40（就是说40%），超过这个比例，整条read就被删除
# -l read小于这个参数设定长度会被丢弃或删除，默认是15,由你的测序bp决定
# -j, --json    输出的json报告文件名，以“.json”结尾
# -h, --html    输出的html报告文件名，以“.html”结尾

4 ###参考基因组下载及建索引（如已有，忽略此步）

## 下载基因组序列
mkdir -p  genome/ARS-UCD1.2  && cd genome/ARS-UCD1.2
nohup wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/002/263/795/GCF_002263795.1_ARS-UCD1.2/GCF_002263795.1_ARS-UCD1.2_genomic.fna.gz &  
gunzip GCF_002263795.1_ARS-UCD1.2_genomic.fna.gz
## 下载gtf文件并解压
wget https://ftp.ncbi.nlm.nih.gov/genomes/refseq/vertebrate_mammalian/Bos_taurus/latest_assembly_versions/GCF_002263795.1_ARS-UCD1.2/GCF_002263795.1_ARS-UCD1.2_genomic.gtf.gz
gunzip GCF_002263795.1_ARS-UCD1.2_genomic.gtf.gz
## 建索引（hisat2软件）
mkdir index & cd index
nohup hisat2-build -p 4 GCF_002263795.1_ARS-UCD1.2_genomic.fna  GCF_002263795.1_ARS-UCD1.2_genomic > hisat2.log &

5 ###reads mapping到参考基因组-----hisat2

mkdir hismap.sam
ls *filter.fastq.gz|cut -d"_" -f 1 |sort -u | while read id ;do (nohup hisat2 -p 8 -x /home/sll/genome-cattle/ARS-UCD1.2/GCF_002263795.1_ARS-UCD1.2_genomic -1 ${id}_1.filter.fastq.gz -2 ${id}_2.filter.fastq.gz -S hismap.sam/${id}.hismap.sam &) ;done


#绵羊参考基因组
/home/sll/genome-sheep/Oar_rambouillet_v1.0-ncbi/GCF_002742125.1_Oar_rambouillet_v1.0_genomic
#牛参考基因组
/home/sll/genome-cattle/ARS-UCD1.2/GCF_002263795.1_ARS-UCD1.2_genomic
# -p: 线程数
# -x: 基因组索引前缀。所下的基因组索引为多个文件，索引前缀到genome为止。
# -1/-2:  fastq输入文件。当输入为单端测序时使用-U 指定输入。单端或双端的输入都可使用逗号分隔输入多个样本，或使用多次-1 -2 / -U 指定多个输入。e.g.: -U sample1.fq.gz,sample2.fq.gz -U sample3.fq.gz
# -S: 输出sam文件路径。

6 ###转为bam文件并排序

cd hisout.sam
ls *sam | cut -d"." -f 1 |sort -u | while read id ;do (nohup samtools view -S ${id}.hismap.sam -b > ${id}.hismap.bam &) ;done
ls *hismap.bam | cut -d"." -f 1 |sort -u | while read id ;do (nohup samtools sort -@ 8 ${id}.hismap.bam -o ${id}_sort.bam &) ;done

# sort: 进行排序
# -@: 线程数
# -o: 输出bam文件
# 最后一项为输入文件

7 ###为sort.bam文件建索引

ls *sort.bam | cut -d"_" -f 1 | sort -u | while read id ;do (nohup samtools index ${id}_sort.bam ${id}_sort.bam.index &) ;done

## index用于建立索引，使之快速访问bam文件

8 ###faturecount计数定量

nohup featureCounts -p -t exon -g gene_id -M -T 8 -a /home/sll/genome-cattle/ARS-UCD1.2/GCF_002263795.1_ARS-UCD1.2_genomic.gtf -o all.featurecounts.txt *_sort.bam &


#绵羊gtf文件
/home/sll/genome-sheep/Oar_rambouillet_v1.0-ncbi/GCF_002742125.1_Oar_rambouillet_v1.0_genomic.gtf
#牛gtf文件
/home/sll/genome-cattle/ARS-UCD1.2/GCF_002263795.1_ARS-UCD1.2_genomic.gtf
#  featureCounts 进行定量，统计比对在这个基因的坐标上的read数
# -t 设置feature-type，-t指定的必须是gtf中有的feature，同时read只有落到这些feature上才会被统计到，默认是“exon”
# -p 双端数据时需要
# -a 转录组注释文件
# -o 输出文件

###报错处理
cat GCF_002263795.1_ARS-UCD1.2_genomic.gtf | grep gene_id | wc -w  #检查原始gtf注释文件gene_id个数
cat GCF_002263795.1_ARS-UCD1.2_genomic.gtf | grep gene_id\ \"\"\; | wc -w  #检查空gene_id值的行数
sed -i -e '/gene_id\ \"\"\;/d' GCF_002263795.1_ARS-UCD1.2_genomic.gtf  #删除包括空gene_id值的行 
cat GCF_002263795.1_ARS-UCD1.2_genomic.gtf | grep gene_id\ \"\"\; | wc -w  #修改后的检查一下 检查包括空gene_id值的行 0为已删除掉

9 ###对featureCounts的结果进行整合，html文件可视化

multiqc all.featurecounts.txt.summary -o  all.counts.summary

10 ###提取定量信息

awk -F '\t' '{print $1,$7,$8,$9,$10,$11,$12}' OFS='\t' all.featurecounts.txt > all_fcount.matrix.txt

# 1列为基因id，7列及以后为样本定量信息
# \t 表示以制表符分割开来

11 ###将矩阵导入R中，采用DESeq2进行差异分析