rnaseq

参考：http://www.bio-info-trainee.com/2218.html

数据来源：https://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP029245&o=acc_s%3Aa

prefetch SRRxxxxx

下载sra
SRAtoolkit也提供了prefetch命令用于下载SRA文件。
格式

ls ../sra/* |while read id; do (fastq-dump --gzip --split-3 -O ./ $id &);done
**sra转fq --split-3**
加上--split-3之后, 会把原来双端拆分成两个文件,但是原来单端并不会保存成两个文件. --gzip就能输出gz格式。

```bash
 ls *gz |xargs fastqc -t 10 -o ../QC/
 ##批量生成质量报告输出到QC文件夹

```bash
 multiqc ./
 ##再QC文件夹下运行multiqc 把所有qc统计为一个

dir='../clean'
cat $1 |while read id
do
	trim_galore -q 25 --phred33 --length 36 -e 0.1 --stringency 3 -o $dir $id
done
##创建一个trim脚本 把fq文件名放入名为1的文件中 

 bash trim.sh 1
 ##运行脚本 把文件1传入

hisat2-build -p 4 hg19.fa hg19
####在hg19参考基因组目录下构建索引

生成8个.ht文件

wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz
tar -zxvf hg19.tar.gz

####或者直接下载索引

(rna) root@PC-jj:/mnt/e/练习/rnaseq/test# ls ../clean/SRR* |while read id ;do(echo $(basename $id)) ; done
#####打印出fq文件名，basename是去掉前面的路径，ls列出的文件名含有路径

(rna) root@PC-jj:/mnt/e/练习/rnaseq/test# ls ../clean/SRR* |while read id ;do(zcat $id |head -1000) > $(basename $id)) ; done
###取前1000行来比对

(rna) root@PC-jj:/mnt/e/练习/rnaseq/test# hisat2 -p 10 -x ../人参考基因组/hg19 -U SRR957677_trimmed.fq.gz  -S align.sam

for i in `seq 77 79`
do 
	hisat2 -t -x ../hg19_index/hg19/genome -U SRR9576${i}_trimmed.fq.gz -S ../aligment/SRR9576${i}.sam &
done
###########对比方法二

---

samtools

---

 samtools view SRR957677.sam
samtools view SRR957677.sam -b >SRR957677.bam

samtools sort SRR957677.bam -o SRR957677_sort.bam &  #按chr排序

samtools index SRR957677_sort.bam &  
#对排序后的序列建立索引，并输出为bai文件，用于快速随机处理。
#在很多情况下，特别是需要显示比对序列的时候，bai文件是必不可少的，
#例如之后的tview命令

head -1000 SRR957677_sort.bam >test_sort.bam

samtools view SRR957677_sort.bam chr1:1234-123456 |less -SN


 samtools flagstat SRR957677_sort.bam
 ##统计bam文件中的比对flag信息，并输出比对统计结果
samtools view -b -f 0x10 SRR957677_sort.bam chr1:1234-123456 >mate_reads.bam

##使用-f或-F参数提取不同匹配情况的read。
#flag是一种描述read比对情况的标记，一种12种，可以搭配使用。
0x1    PAIRED    paired-end (or multiple-segment) sequencing technology
0x2    PROPER_PAIR         each segment properly aligned according to the aligner
0x4    UNMAP               segment unmapped
0x8    MUNMAP              next segment in the template unmapped
0x10   REVERSE             SEQ is reverse complemented
0x20   MREVERSE            SEQ of the next segment in the template is reverse complemented
0x40   READ1               the first segment in the template
0x80   READ2               the last segment in the template
0x100  SECONDARY           secondary alignment
0x200  QCFAIL              not passing quality controls
0x400  DUP                 PCR or optical duplicate
0x800  SUPPLEMENTARY       supplementary alignment


想用samtools筛选恰好配对的read,就需要用0x10

6. tview

作用
以文本定位查看器的方式来展示各条reads的比对情况
格式
$ samtools tview [options] <排序后的bam文件> [参考基因组fasta文件]

htseq

htseq-count -f bam SRR957678_sort.bam  ../homo/gencode.v19.annotation.gtf.gz >SRR957678.count

###计数

 wc -l ls *.count
paste *.count |cut -f 1,2,4,6,8 --output-delimiter $'\t' >countmatrix
##--output-delimiter $'\t'  指定输出分隔符为\t

 cat countmatrix |cut -f 1 |cut -d . -f 1 |head
 ##去除geneid名称的小数点，
 #gencode的注释文件中的gene_id（如ENSG00000000003.10）
#在EBI是不能搜索到的，只保留ENSG00000000003这部分

 paste gene_id countmatrix |cut -f 1,3,4,5,6 --output-delimiter $'\t' > ../DEG/countmatrix1
 ###替换原来含小数点的geneid